Repository: PennyLaneAI/pennylane
Pull request #1102 (2021-02-22 16:26:36)

Add kernel module
This PR adds a new submodule to Pennylane, `qml.kernels`. It contains the following functionalities:
* Create Quantum Embedding Kernels from an ansatz template via the `EmbeddingKernel` class
* Convenience methods for the calculation of kernel matrices, geared towards use with scikit-learn
* Post-processing to mitigate non-positive semidefiniteness of kernel matrices
* Quality measures for kernels relative to datasets, namely polarization and kernel-target alignment

Additionally, tests for the relevant new functionalities were added.

codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102?src=pr&el=h1) Report
> Merging [#1102](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102?src=pr&el=desc) (3c89b70) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/c8673eb45d9b3d435867f04b9943a5c4b53322be?el=desc) (c8673eb) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master    #1102   +/-   ##
=======================================
  Coverage   97.73%   97.73%           
=======================================
  Files         154      155    +1     
  Lines       11681    11683    +2     
=======================================
+ Hits        11416    11418    +2     
  Misses        265      265           
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pennylane/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102/diff?src=pr&el=tree#diff-cGVubnlsYW5lL19faW5pdF9fLnB5) | `97.14% <100.00%> (+0.04%)` | :arrow_up: |
| [pennylane/kernels/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102/diff?src=pr&el=tree#diff-cGVubnlsYW5lL2tlcm5lbHMvX19pbml0X18ucHk=) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102?src=pr&el=footer). Last update [c8673eb...3c89b70](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1102?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

johannesjmeyer: Call results:
- Don't fix signatures
- Have utilities in kernels module
- Remove Kernel class
- cvxpy as a "when-needed" import, optional dependency as in pip install pennylane [kernels]
johannesjmeyer: @mariaschuld @josh146 @co9olguy @antalszava, you can have a first go at the review, maybe first see if you are happy with the structure. @dwierichs and me sat down to implement the changes as per our meeting last week. We moved one method to utils, removed the object-oriented interface, import cvxpy on a "when needed" basis, added a "kernels" setup option that installs cvxpy and improved the docstrings. 

Usage can be checked out in the corresponding demo at https://github.com/PennyLaneAI/qml/pull/266. 
dwierichs: Thanks for the detailed feedback! :)

I implemented all changes except for 
1. moving the literature links to the docstrings (waiting for Johannes' response)
2. changing the starting year of the copyright period (i.e. I went with 2018-2021 for now)

Regarding the type hint, I did not change anything, is that alright?
@mariaschuld, maybe you can have a look at my answer to your comment on using symmetry in `kernel_matrix` above :)
mariaschuld: @dwierichs thanks for addressing the comments; let me know when everything is passing and ready for final review, should be easy to approve then!
dwierichs: Looks like this could be merged :) @mariaschuld @josh146 
josh146: @dwierichs looks like there are a couple of rendering issues in the [new doc page](https://pennylane--1102.org.readthedocs.build/en/1102/code/qml_kernels.html)

Note: you can access the built docs for the PR by clicking on 'details' next to the 'readthedocs' CI check:

![image](https://user-images.githubusercontent.com/2959003/118661790-078cd100-b822-11eb-9c50-91a3f7563bc5.png)

josh146: Oh, one more thing -- apologies for just catching this. All PL docstrings must have example usage, to provide a minimal example of how the function is used. Check out our PL documentation guidelines here: https://pennylane.readthedocs.io/en/stable/development/guide/documentation.html#code-examples
johannesjmeyer: I added some final touches, the only big one being to move `kernel_matrix` and `square_kernel_matrix` to a new `kernels.utils` module because they didn't really fit into `cost_functions` and to rename `kernel_polarization` to the correct `kernel_polarity`.

Otherwise this is probably very close to merging, thanks a ton to @dwierichs who did all the heavy lifting for getting this shipshape!
dwierichs: Is there something to be done by us for the codecov to go up? It seems to be cause by latest merges into master?

johannesjmeyer: @josh146, will you push the button? :) 


Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 25ce7e5a0..8710e9349 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -1,6 +1,49 @@
 # Release 0.16.0-dev (development release)
 
 <h3>New features since last release</h3>
+
+* PennyLane now has a ``kernels`` module.
+  It provides basic functionalities for working with quantum kernels as well as 
+  post-processing methods to mitigate sampling errors and device noise:
+
+```python
+import pennylane as qml
+from pennylane import numpy as np
+
+num_wires = 6
+wires = range(num_wires)
+
+dev = qml.device('default.qubit', wires=num_wires)
+
+@qml.qnode(dev)
+def kernel_circuit(x1, x2):
+    qml.templates.AngleEmbedding(x1, wires=wires)
+    qml.adjoint(qml.templates.AngleEmbedding)(x2, wires=wires)
+    return qml.probs(wires)
+
+kernel = lambda x1, x2: kernel_circuit(x1, x2)[0]
+
+# "Training feature vectors"
+X_train = np.random.random((10, 6))
+# Create symmetric square kernel matrix (for training)
+K = qml.kernels.square_kernel_matrix(X_train, kernel)
+# Add some (symmetric) Gaussian noise to the kernel matrix.
+N = np.random.randn(10, 10)
+K += (N + N.T) / 2
+        
+K1 = qml.kernels.displace_matrix(K)
+K2 = qml.kernels.closest_psd_matrix(K)
+K3 = qml.kernels.threshold_matrix(K)
+K4 = qml.kernels.mitigate_depolarizing_noise(K, num_wires, method='single')
+K5 = qml.kernels.mitigate_depolarizing_noise(K, num_wires, method='average')
+K6 = qml.kernels.mitigate_depolarizing_noise(K, num_wires, method='split_channel')
+
+# "Testing feature vectors"
+X_test = np.random.random((5, 6))
+# Compute kernel between test and training data.
+K_test = qml.kernels.kernel_matrix(X_train, X_test, kernel)
+```
+
 * Added CPhase operation as an alias for ControlledPhaseShift operation
   [(#1319)](https://github.com/PennyLaneAI/pennylane/pull/1319).
 
@@ -307,8 +350,8 @@ random_mat2 = rng.standard_normal(3, requires_grad=False)
 This release contains contributions from (in alphabetical order):
 
 Marius Aglitoiu, Vishnu Ajith, Thomas Bromley, Jack Ceroni, Miruna Daian, Olivia Di Matteo,
-Tanya Garg, Diego Guala, Anthony Hayes, Josh Izaac, Pavan Jayasinha, Ryan Levy, Nahum Sá,
-Brian Shi, Antal Száva, Vincent Wong, Alberto Maldonado.
+Tanya Garg, Diego Guala, Anthony Hayes, Josh Izaac, Pavan Jayasinha, Ryan Levy, Nahum Sá, Johannes Jakob Meyer,
+Brian Shi, Antal Száva, David Wierichs, Vincent Wong, Alberto Maldonado.
 
 # Release 0.15.1 (current release)
 
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 25a8e4011..0a6883ed6 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -11,5 +11,5 @@ jobs:
       with:
         docs-folder: "doc/"
         repo-token: "${{ secrets.GITHUB_TOKEN }}"
-        pre-build-command: "apt install -y graphviz && pip3 install -r requirements.txt && pip3 install . && pip3 install -r qchem/requirements.txt && pip3 install ./qchem"
+        pre-build-command: "apt install -y graphviz && pip3 install -r doc/requirements.txt && pip3 install . && pip3 install -r qchem/requirements.txt && pip3 install ./qchem"
         build-command: "sphinx-build -b html . _build -W --keep-going"
diff --git a/doc/code/qml_kernels.rst b/doc/code/qml_kernels.rst
new file mode 100644
index 000000000..43ee81091
--- /dev/null
+++ b/doc/code/qml_kernels.rst
@@ -0,0 +1,9 @@
+qml.kernels
+===========
+
+.. currentmodule:: pennylane.kernels
+
+.. automodapi:: pennylane.kernels
+    :no-heading:
+    :include-all-objects:
+    :no-inheritance-diagram:
diff --git a/doc/development/guide/documentation.rst b/doc/development/guide/documentation.rst
index 43156d0d5..4cc3127a9 100644
--- a/doc/development/guide/documentation.rst
+++ b/doc/development/guide/documentation.rst
@@ -70,7 +70,7 @@ Punctuation and spelling
           ``QNodeCollection`` should be executed. This can either be a single device, or a list
           of devices of length ``len(observables)``.
 
-* Comments may written more informally then docstrings, as long as consistency and clarity is maintained.
+* Comments may be written more informally than docstrings, as long as consistency and clarity are maintained.
   Capitalization and punctuation should be used with multi-sentence comments to aid with readability.
 
 
diff --git a/doc/index.rst b/doc/index.rst
index af99c30c5..dfee13bfc 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -207,6 +207,7 @@ PennyLane is **free** and **open source**, released under the Apache License, Ve
    code/qml_operation
    code/qml_devices
    code/qml_grouping
+   code/qml_kernels
    code/qml_math
    code/qml_qaoa
    code/qml_qchem
diff --git a/doc/requirements.txt b/doc/requirements.txt
index 4e0399704..a3fe65cc5 100644
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@@ -2,6 +2,7 @@
 pip
 appdirs
 autograd
+autoray
 jax==0.2.7
 jaxlib==0.1.57
 m2r2
diff --git a/pennylane/__init__.py b/pennylane/__init__.py
index 168ecc1ad..3b7f8adc5 100644
--- a/pennylane/__init__.py
+++ b/pennylane/__init__.py
@@ -22,6 +22,7 @@ import pkg_resources
 from semantic_version import Spec, Version
 
 import pennylane.init
+import pennylane.kernels
 import pennylane.math
 import pennylane.operation
 import pennylane.qaoa as qaoa
diff --git a/pennylane/kernels/__init__.py b/pennylane/kernels/__init__.py
new file mode 100644
index 000000000..5f10e05cf
--- /dev/null
+++ b/pennylane/kernels/__init__.py
@@ -0,0 +1,82 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+r"""
+This subpackage defines functions that relate to quantum kernel methods.
+On one hand this includes functions to call a quantum kernel systematically
+on training and test datasets to obtain the *kernel matrix*.
+On the other hand it provides postprocessing methods for those kernel
+matrices which can be used to mitigate device noise and sampling errors.
+
+Given a kernel
+
+.. math ::
+
+    k: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}, \quad
+    (x_1, x_2)\mapsto k(x_1, x_2)
+
+the kernel matrix of :math:`k` on a training dataset
+:math:`\{(x_1, y_1),\cdots (x_n, y_n)\}` with :math:`x_i\in\mathbb{R}^d`
+and :math:`y_i\in\{-1, 1\}` is defined as
+
+.. math ::
+
+    K_{ij} = k(x_i, x_j).
+
+For valid kernels, this is a real symmetric positive semi-definite matrix.
+We also define the *ideal kernel matrix* for the training dataset which
+perfectly predicts whether two points have identical labels or not:
+
+.. math ::
+
+    K^\ast_{ij} = y_i y_j
+
+We can measure the similarity between :math:`K` and :math:`K^\ast`,
+through the *kernel polarity* which can be expressed as the Frobenius inner
+product between the two matrices:
+
+.. math ::
+
+    \operatorname{P}(k) = \langle K^\ast, K \rangle_F = \sum_{i,j=1}^n y_i y_j k(x_i, x_j)
+
+Additionally, there is the *kernel-target alignment*, which is the normalized
+counterpart to the kernel polarity:
+
+.. math ::
+
+    \operatorname{TA}(k) &= \frac{P(k)}{\lVert K^\ast \rVert_F\;\lVert K \rVert_F}\\
+        \lVert K\rVert_F &= \sqrt{\sum_{i,j=1}^n k(x_i, x_j)^2}\\
+   \lVert K^\ast\rVert_F &= \sqrt{\sum_{i,j=1}^n (y_iy_j)^2}
+
+For datasets with different numbers of training points per class the labels are rescaled
+by the number of datapoints in the respective class to avoid that kernel polarity and
+kernel-target alignment are dominated by the properties of the kernel for just a single class.
+
+Given a callable kernel function, all these quantities can readily be computed
+using the methods in this module.
+"""
+from .cost_functions import (
+    polarity,
+    target_alignment,
+)
+from .postprocessing import (
+    threshold_matrix,
+    displace_matrix,
+    flip_matrix,
+    closest_psd_matrix,
+    mitigate_depolarizing_noise,
+)
+from .utils import (
+    kernel_matrix,
+    square_kernel_matrix,
+)
diff --git a/pennylane/kernels/cost_functions.py b/pennylane/kernels/cost_functions.py
new file mode 100644
index 000000000..b4abd8e04
--- /dev/null
+++ b/pennylane/kernels/cost_functions.py
@@ -0,0 +1,177 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This file contains functionalities for kernel related costs.
+See `here <https://www.doi.org/10.1007/s10462-012-9369-4>`_ for a review.
+"""
+from pennylane import numpy as np
+from ..utils import frobenius_inner_product
+from .utils import square_kernel_matrix
+
+
+def polarity(
+    X,
+    Y,
+    kernel,
+    assume_normalized_kernel=False,
+    rescale_class_labels=True,
+    normalize=False,
+):
+    r"""Polarity of a given kernel function.
+
+    For a dataset with feature vectors :math:`\{x_i\}` and associated labels :math:`\{y_i\}`,
+    the polarity of the kernel function :math:`k` is given by
+
+    .. math ::
+
+        \operatorname{P}(k) = \sum_{i,j=1}^n y_i y_j k(x_i, x_j)
+
+    If the dataset is unbalanced, that is if the numbers of datapoints in the
+    two classes :math:`n_+` and :math:`n_-` differ,
+    ``rescale_class_labels=True`` will apply a rescaling according to
+    :math:`\tilde{y}_i = \frac{y_i}{n_{y_i}}`. This is activated by default
+    and only results in a prefactor that depends on the size of the dataset
+    for balanced datasets.
+
+    The keyword argument ``assume_normalized_kernel`` is passed to
+    :func:`~.kernels.square_kernel_matrix`, for the computation
+    :func:`~.utils.frobenius_inner_product` is used.
+
+    Args:
+        X (list[datapoint]): List of datapoints.
+        Y (list[float]): List of class labels of datapoints, assumed to be either -1 or 1.
+        kernel ((datapoint, datapoint) -> float): Kernel function that maps datapoints to kernel value.
+        assume_normalized_kernel (bool, optional): Assume that the kernel is normalized, i.e.
+            the kernel evaluates to 1 when both arguments are the same datapoint.
+        rescale_class_labels (bool, optional): Rescale the class labels. This is important to take
+            care of unbalanced datasets.
+        normalize (bool): If True, rescale the polarity to the target_alignment.
+
+    Returns:
+        float: The kernel polarity.
+
+    **Example:**
+
+    Consider a simple kernel function based on :class:`~.templates.embeddings.AngleEmbedding`:
+
+    .. code-block :: python
+
+        dev = qml.device('default.qubit', wires=2, shots=None)
+        @qml.qnode(dev)
+        def circuit(x1, x2):
+            qml.templates.AngleEmbedding(x1, wires=dev.wires)
+            qml.adjoint(qml.templates.AngleEmbedding)(x2, wires=dev.wires)
+            return qml.probs(wires=dev.wires)
+
+        kernel = lambda x1, x2: circuit(x1, x2)[0]
+
+    We can then compute the polarity on a set of 4 (random) feature
+    vectors ``X`` with labels ``Y`` via
+
+    >>> X = np.random.random((4, 2))
+    >>> Y = np.array([-1, -1, 1, 1])
+    >>> qml.kernels.polarity(X, Y, kernel)
+    tensor(0.04361349, requires_grad=True)
+    """
+    K = square_kernel_matrix(X, kernel, assume_normalized_kernel=assume_normalized_kernel)
+
+    if rescale_class_labels:
+        nplus = np.count_nonzero(np.array(Y) == 1)
+        nminus = len(Y) - nplus
+        _Y = np.array([y / nplus if y == 1 else y / nminus for y in Y])
+    else:
+        _Y = np.array(Y)
+
+    T = np.outer(_Y, _Y)
+
+    return frobenius_inner_product(K, T, normalize=normalize)
+
+
+def target_alignment(
+    X,
+    Y,
+    kernel,
+    assume_normalized_kernel=False,
+    rescale_class_labels=True,
+):
+    r"""Target alignment of a given kernel function.
+
+    This function is an alias for :func:`~.kernels.polarity` with ``normalize=True``.
+
+    For a dataset with feature vectors :math:`\{x_i\}` and associated labels :math:`\{y_i\}`, the
+    target alignment of the kernel function :math:`k` is given by
+
+    .. math ::
+
+        \operatorname{TA}(k) = \frac{\sum_{i,j=1}^n y_i y_j k(x_i, x_j)}
+        {\sqrt{\sum_{i,j=1}^n y_i y_j} \sqrt{\sum_{i,j=1}^n k(x_i, x_j)^2}}
+
+    If the dataset is unbalanced, that is if the numbers of datapoints in the
+    two classes :math:`n_+` and :math:`n_-` differ,
+    ``rescale_class_labels=True`` will apply a rescaling according to
+    :math:`\tilde{y}_i = \frac{y_i}{n_{y_i}}`. This is activated by default
+    and only results in a prefactor that depends on the size of the dataset
+    for balanced datasets.
+
+    Args:
+        X (list[datapoint]): List of datapoints
+        Y (list[float]): List of class labels of datapoints, assumed to be either -1 or 1.
+        kernel ((datapoint, datapoint) -> float): Kernel function that maps datapoints to kernel value.
+        assume_normalized_kernel (bool, optional): Assume that the kernel is normalized, i.e.
+            the kernel evaluates to 1 when both arguments are the same datapoint.
+        rescale_class_labels (bool, optional): Rescale the class labels. This is important to take
+            care of unbalanced datasets.
+
+    Returns:
+        float: The kernel-target alignment.
+
+    **Example:**
+
+    Consider a simple kernel function based on :class:`~.templates.embeddings.AngleEmbedding`:
+
+    .. code-block :: python
+
+        dev = qml.device('default.qubit', wires=2, shots=None)
+        @qml.qnode(dev)
+        def circuit(x1, x2):
+            qml.templates.AngleEmbedding(x1, wires=dev.wires)
+            qml.adjoint(qml.templates.AngleEmbedding)(x2, wires=dev.wires)
+            return qml.probs(wires=dev.wires)
+
+        kernel = lambda x1, x2: circuit(x1, x2)[0]
+
+    We can then compute the kernel-target alignment on a set of 4 (random)
+    feature vectors ``X`` with labels ``Y`` via
+
+    >>> X = np.random.random((4, 2))
+    >>> Y = np.array([-1, -1, 1, 1])
+    >>> qml.kernels.target_alignment(X, Y, kernel)
+    tensor(0.01124802, requires_grad=True)
+
+    We can see that this is equivalent to using ``normalize=True`` in
+    ``polarity``:
+
+    >>> target_alignment = qml.kernels.target_alignment(X, Y, kernel)
+    >>> normalized_polarity = qml.kernels.polarity(X, Y, kernel, normalize=True)
+    >>> np.isclose(target_alignment, normalized_polarity)
+    tensor(True, requires_grad=True)
+    """
+    return polarity(
+        X,
+        Y,
+        kernel,
+        assume_normalized_kernel=assume_normalized_kernel,
+        rescale_class_labels=rescale_class_labels,
+        normalize=True,
+    )
diff --git a/pennylane/kernels/postprocessing.py b/pennylane/kernels/postprocessing.py
new file mode 100644
index 000000000..70b5aec01
--- /dev/null
+++ b/pennylane/kernels/postprocessing.py
@@ -0,0 +1,300 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This file contains functionalities for postprocessing of kernel matrices.
+"""
+from pennylane import numpy as np
+
+
+def threshold_matrix(K):
+    r"""Remove negative eigenvalues from the given kernel matrix.
+
+    This method yields the closest positive semi-definite matrix in
+    any unitarily invariant norm, e.g. the Frobenius norm.
+
+    Args:
+        K (array[float]): Kernel matrix, assumed to be symmetric.
+
+    Returns:
+        array[float]: Kernel matrix with cropped negative eigenvalues.
+
+    **Example:**
+
+    Consider a symmetric matrix with both positive and negative eigenvalues:
+
+    .. code-block :: pycon
+
+        >>> K = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 2]])
+        >>> np.linalg.eigvalsh(K)
+        array([-1.,  1.,  2.])
+
+    We then can threshold/truncate the eigenvalues of the matrix via
+
+    .. code-block :: pycon
+
+        >>> K_thresh = qml.kernels.threshold_matrix(K)
+        >>> np.linalg.eigvalsh(K_thresh)
+        array([0.,  1.,  2.])
+
+    If the input matrix does not have negative eigenvalues, ``threshold_matrix``
+    does not have any effect.
+    """
+    w, v = np.linalg.eigh(K)
+
+    if w[0] < 0:
+        # Transform spectrum: Threshold/clip at 0.
+        w0 = np.clip(w, 0, None)
+
+        return (v * w0) @ np.transpose(v)
+
+    return K
+
+
+def displace_matrix(K):
+    r"""Remove negative eigenvalues from the given kernel matrix by adding a multiple
+    of the identity matrix.
+
+    This method keeps the eigenvectors of the matrix intact.
+
+    Args:
+        K (array[float]): Kernel matrix, assumed to be symmetric.
+
+    Returns:
+        array[float]: Kernel matrix with eigenvalues offset by adding the identity.
+
+    **Example:**
+
+    Consider a symmetric matrix with both positive and negative eigenvalues:
+
+    .. code-block :: pycon
+
+        >>> K = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 2]])
+        >>> np.linalg.eigvalsh(K)
+        array([-1.,  1.,  2.])
+
+    We then can shift all eigenvalues of the matrix by adding the identity matrix
+    multiplied with the absolute value of the smallest (the most negative, that is)
+    eigenvalue:
+
+    .. code-block :: pycon
+
+        >>> K_displaced = qml.kernels.displace_matrix(K)
+        >>> np.linalg.eigvalsh(K_displaced)
+        array([0.,  2.,  3.])
+
+    If the input matrix does not have negative eigenvalues, ``displace_matrix``
+    does not have any effect.
+    """
+    wmin = np.linalg.eigvalsh(K)[0]
+
+    if wmin < 0:
+        return K - np.eye(K.shape[0]) * wmin
+
+    return K
+
+
+def flip_matrix(K):
+    r"""Remove negative eigenvalues from the given kernel matrix by taking the absolute value.
+
+    This method keeps the eigenvectors of the matrix intact.
+
+    Args:
+        K (array[float]): Kernel matrix, assumed to be symmetric.
+
+    Returns:
+        array[float]: Kernel matrix with flipped negative eigenvalues.
+
+    Reference:
+        This method is introduced in `arXiv:2103.16774 <https://arxiv.org/abs/2103.16774>`_.
+
+    **Example:**
+
+    Consider a symmetric matrix with both positive and negative eigenvalues:
+
+    .. code-block :: pycon
+
+        >>> K = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 2]])
+        >>> np.linalg.eigvalsh(K)
+        array([-1.,  1.,  2.])
+
+    We then can invert the sign of all negative eigenvalues of the matrix, obtaining
+    non-negative eigenvalues only:
+
+    .. code-block :: pycon
+
+        >>> K_flipped = qml.kernels.flip_matrix(K)
+        >>> np.linalg.eigvalsh(K_flipped)
+        array([1.,  1.,  2.])
+
+    If the input matrix does not have negative eigenvalues, ``flip_matrix``
+    does not have any effect.
+    """
+    w, v = np.linalg.eigh(K)
+
+    if w[0] < 0:
+        # Transform spectrum: absolute value
+        w_abs = np.abs(w)
+
+        return (v * w_abs) @ np.transpose(v)
+
+    return K
+
+
+def closest_psd_matrix(K, fix_diagonal=False, solver=None, **kwargs):
+    r"""Return the closest positive semi-definite matrix to the given kernel matrix.
+
+    This method either fixes the diagonal entries to be 1
+    (``fix_diagonal=True``) or keeps the eigenvectors intact (``fix_diagonal=False``),
+    in which case it reduces to the method :func:`~.kernels.threshold_matrix`.
+    For ``fix_diagonal=True`` a semi-definite program is solved.
+
+    Args:
+        K (array[float]): Kernel matrix, assumed to be symmetric.
+        fix_diagonal (bool): Whether to fix the diagonal to 1.
+        solver (str, optional): Solver to be used by cvxpy. Defaults to CVXOPT.
+        kwargs (kwarg dict): Passed to cvxpy.Problem.solve().
+
+    Returns:
+        array[float]: closest positive semi-definite matrix in Frobenius norm.
+
+    Comments:
+        Requires cvxpy and the used solver (default CVXOPT) to be installed if ``fix_diagonal=True``.
+
+    Reference:
+        This method is introduced in `arXiv:2105.02276 <https://arxiv.org/abs/2105.02276>`_.
+
+    **Example:**
+
+    Consider a symmetric matrix with both positive and negative eigenvalues:
+
+    .. code-block :: pycon
+
+        >>> K = np.array([[0.9, 1.], [1., 0.9]])
+        >>> np.linalg.eigvalsh(K)
+        array([-0.1, 1.9])
+
+    The positive semi-definite matrix that is closest to this matrix in any unitarily
+    invariant norm is then given by the matrix with the eigenvalues thresholded at 0,
+    as computed by :func:`~.kernels.threshold_matrix`:
+
+    .. code-block :: pycon
+
+        >>> K_psd = qml.kernels.closest_psd_matrix(K)
+        >>> K_psd
+        tensor([[0.95, 0.95],
+                [0.95, 0.95]], requires_grad=True)
+        >>> np.linalg.eigvalsh(K_psd)
+        array([0., 1.9])
+        >>> np.allclose(K_psd, qml.kernels.threshold_matrix(K))
+        True
+
+    However, for quantum kernel matrices we may want to restore the value 1 on the
+    diagonal:
+
+    .. code-block :: pycon
+
+        >>> K_psd = qml.kernels.closest_psd_matrix(K, fix_diagonal=True)
+        >>> K_psd
+        array([[1.        , 0.99998008],
+               [0.99998008, 1.        ]])
+        >>> np.linalg.eigvalsh(K_psd)
+        array([1.99162415e-05, 1.99998008e+00])
+
+    If the input matrix does not have negative eigenvalues and ``fix_diagonal=False``,
+    ``closest_psd_matrix`` does not have any effect.
+    """
+    if not fix_diagonal:
+        return threshold_matrix(K)
+    try:
+        import cvxpy as cp
+
+        if solver is None:
+            solver = cp.CVXOPT
+    except ImportError as e:
+        raise ImportError("CVXPY is required for this post-processing method.") from e
+
+    X = cp.Variable(K.shape, PSD=True)
+    constraint = [cp.diag(X) == 1.0] if fix_diagonal else []
+    objective_fn = cp.norm(X - K, "fro")
+    problem = cp.Problem(cp.Minimize(objective_fn), constraint)
+
+    try:
+        problem.solve(solver=solver, **kwargs)
+    except Exception:
+        try:
+            problem.solve(solver=solver, verbose=True, **kwargs)
+        except Exception as e:
+            raise RuntimeError("CVXPY solver did not converge.") from e
+
+    return X.value
+
+
+def mitigate_depolarizing_noise(K, num_wires, method, use_entries=None):
+    r"""Estimate depolarizing noise rate(s) using on the diagonal entries of a kernel
+    matrix and mitigate the noise, assuming a global depolarizing noise model.
+
+    Args:
+        K (array[float]): Noisy kernel matrix.
+        num_wires (int): Number of wires/qubits of the quantum embedding kernel.
+        method (``'single'`` | ``'average'`` | ``'split_channel'``): Strategy for mitigation
+
+            * ``'single'``: An alias for ``'average'`` with ``len(use_entries)=1``.
+            * ``'average'``: Estimate a global noise rate based on the average of the diagonal
+              entries in ``use_entries``, which need to be measured on the quantum computer.
+            * ``'split_channel'``: Estimate individual noise rates per embedding, requiring
+              all diagonal entries to be measured on the quantum computer.
+        use_entries (array[int]): Diagonal entries to use if method in ``['single', 'average']``.
+            If ``None``, defaults to ``[0]`` (``'single'``) or ``range(len(K))`` (``'average'``).
+
+    Returns:
+        array[float]: Mitigated kernel matrix.
+
+    Reference:
+        This method is introduced in Section V in
+        `arXiv:2105.02276 <https://arxiv.org/abs/2105.02276>`_.
+
+    **Example:**
+
+    For an example usage of ``mitigate_depolarizing_noise`` please refer to the
+    `PennyLane demo on the kernel module <https://github.com/PennyLaneAI/qml/tree/master/demonstrations/tutorial_kernel_module.py>`_ or `the postprocessing demo for arXiv:2105.02276 <https://github.com/thubregtsen/qhack/blob/master/paper/post_processing_demo.py>`_.
+    """
+    dim = 2 ** num_wires
+
+    if method == "single":
+        if use_entries is None:
+            use_entries = (0,)
+        diagonal_element = K[use_entries[0], use_entries[0]]
+        noise_rate = (1 - diagonal_element) * dim / (dim - 1)
+        mitigated_matrix = (K - noise_rate / dim) / (1 - noise_rate)
+
+    elif method == "average":
+        if use_entries is None:
+            diagonal_elements = np.diag(K)
+        else:
+            diagonal_elements = np.diag(K)[np.array(use_entries)]
+        noise_rates = (1 - diagonal_elements) * dim / (dim - 1)
+        mean_noise_rate = np.mean(noise_rates)
+        mitigated_matrix = (K - mean_noise_rate / dim) / (1 - mean_noise_rate)
+
+    elif method == "split_channel":
+        eff_noise_rates = np.clip((1 - np.diag(K)) * dim / (dim - 1), 0.0, 1.0)
+        noise_rates = 1 - np.sqrt(1 - eff_noise_rates)
+        inverse_noise = (
+            -np.outer(noise_rates, noise_rates)
+            + noise_rates.reshape((1, len(K)))
+            + noise_rates.reshape((len(K), 1))
+        )
+        mitigated_matrix = (K - inverse_noise / dim) / (1 - inverse_noise)
+
+    return mitigated_matrix
diff --git a/pennylane/kernels/utils.py b/pennylane/kernels/utils.py
new file mode 100644
index 000000000..89535445d
--- /dev/null
+++ b/pennylane/kernels/utils.py
@@ -0,0 +1,122 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This file contains functionalities that simplify working with kernels.
+"""
+from pennylane import numpy as np
+
+
+def square_kernel_matrix(X, kernel, assume_normalized_kernel=False):
+    r"""Computes the square matrix of pairwise kernel values for a given dataset.
+
+    Args:
+        X (list[datapoint]): List of datapoints
+        kernel ((datapoint, datapoint) -> float): Kernel function that maps
+            datapoints to kernel value.
+        assume_normalized_kernel (bool, optional): Assume that the kernel is normalized, in
+            which case the diagonal of the kernel matrix is set to 1, avoiding unnecessary
+            computations.
+
+    Returns:
+        array[float]: The square matrix of kernel values.
+
+    **Example:**
+
+    Consider a simple kernel function based on :class:`~.templates.embeddings.AngleEmbedding`:
+
+    .. code-block :: python
+
+        dev = qml.device('default.qubit', wires=2, shots=None)
+        @qml.qnode(dev)
+        def circuit(x1, x2):
+            qml.templates.AngleEmbedding(x1, wires=dev.wires)
+            qml.adjoint(qml.templates.AngleEmbedding)(x2, wires=dev.wires)
+            return qml.probs(wires=dev.wires)
+
+        kernel = lambda x1, x2: circuit(x1, x2)[0]
+
+    We can then compute the kernel matrix on a set of 4 (random) feature
+    vectors ``X`` via
+
+    >>> X = np.random.random((4, 2))
+    >>> qml.kernels.square_kernel_matrix(X, kernel)
+    tensor([[1.        , 0.9532702 , 0.96864001, 0.90932897],
+            [0.9532702 , 1.        , 0.99727485, 0.95685561],
+            [0.96864001, 0.99727485, 1.        , 0.96605621],
+            [0.90932897, 0.95685561, 0.96605621, 1.        ]], requires_grad=True)
+    """
+    N = len(X)
+    matrix = [0] * N ** 2
+
+    for i in range(N):
+        for j in range(i, N):
+            if assume_normalized_kernel and i == j:
+                matrix[N * i + j] = 1.0
+            else:
+                matrix[N * i + j] = kernel(X[i], X[j])
+                matrix[N * j + i] = matrix[N * i + j]
+
+    return np.array(matrix).reshape((N, N))
+
+
+def kernel_matrix(X1, X2, kernel):
+    r"""Computes the matrix of pairwise kernel values for two given datasets.
+
+    Args:
+        X1 (list[datapoint]): List of datapoints (first argument)
+        X2 (list[datapoint]): List of datapoints (second argument)
+        kernel ((datapoint, datapoint) -> float): Kernel function that maps datapoints to kernel value.
+
+    Returns:
+        array[float]: The square matrix of kernel values.
+
+    **Example:**
+
+    Consider a simple kernel function based on :class:`~.templates.embeddings.AngleEmbedding`:
+
+    .. code-block :: python
+
+        dev = qml.device('default.qubit', wires=2, shots=None)
+        @qml.qnode(dev)
+        def circuit(x1, x2):
+            qml.templates.AngleEmbedding(x1, wires=dev.wires)
+            qml.adjoint(qml.templates.AngleEmbedding)(x2, wires=dev.wires)
+            return qml.probs(wires=dev.wires)
+
+        kernel = lambda x1, x2: circuit(x1, x2)[0]
+
+    With this method we can systematically evaluate the kernel function ``kernel`` on
+    pairs of datapoints, where the points stem from different datasets, like a training
+    and a test dataset.
+
+    >>> X_train = np.random.random((4,2))
+    >>> X_test = np.random.random((3,2))
+    >>> qml.kernels.kernel_matrix(X_train, X_test, kernel)
+    tensor([[0.88875298, 0.90655175, 0.89926447],
+            [0.93762197, 0.98163781, 0.93076383],
+            [0.91977339, 0.9799841 , 0.91582698],
+            [0.80376818, 0.98720925, 0.79349212]], requires_grad=True)
+
+    As we can see, for :math:`n` and :math:`m` datapoints in the first and second
+    dataset respectively, the output matrix has the shape :math:`n\times m`.
+    """
+    N = len(X1)
+    M = len(X2)
+
+    matrix = [0] * N * M
+    for i in range(N):
+        for j in range(M):
+            matrix[M * i + j] = kernel(X1[i], X2[j])
+
+    return np.array(matrix).reshape((N, M))
diff --git a/pennylane/utils.py b/pennylane/utils.py
index 2cbe58ef4..8ee746eaf 100644
--- a/pennylane/utils.py
+++ b/pennylane/utils.py
@@ -454,3 +454,37 @@ def expand_vector(vector, original_wires, expanded_wires):
     expanded_tensor = np.moveaxis(expanded_tensor, original_indices, wire_indices)
 
     return expanded_tensor.reshape(2 ** M)
+
+
+def frobenius_inner_product(A, B, normalize=False):
+    r"""Frobenius inner product between two matrices.
+
+    .. math::
+
+        \langle A, B \rangle_F = \sum_{i,j=1}^n A_{ij} B_{ij} = \operatorname{tr} (A^T B)
+
+    The Frobenius inner product is equivalent to the Hilbert-Schmidt inner product for
+    matrices with real-valued entries.
+
+    Args:
+        A (array[float]): First matrix, assumed to be a square array.
+        B (array[float]): Second matrix, assumed to be a square array.
+        normalize (bool): If True, divide the inner_product by the Frobenius norms of A and B.
+            Defaults to False.
+
+    Returns:
+        float: Frobenius inner product of A and B
+
+    **Example**
+
+    >>> A = np.random.random((3,3))
+    >>> B = np.random.random((3,3))
+    >>> qml.utils.frobenius_inner_product(A, B)
+    3.091948202943376
+    """
+    inner_product = np.sum(A * B)
+
+    if normalize:
+        inner_product /= np.linalg.norm(A, "fro") * np.linalg.norm(B, "fro")
+
+    return inner_product
diff --git a/requirements.txt b/requirements.txt
index 04e6b254c..182e935d1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,7 @@
 numpy
 scipy
+cvxpy
+cvxopt
 networkx
 tensornetwork==0.3
 autograd
diff --git a/setup.py b/setup.py
index 4e414deb0..5531f525d 100644
--- a/setup.py
+++ b/setup.py
@@ -58,6 +58,7 @@ info = {
     'long_description_content_type': "text/markdown",
     'provides': ["pennylane"],
     'install_requires': requirements,
+    'extras_require': {'kernels': ['cvxpy', 'cvxopt']},
     'package_data': {'pennylane': ['devices/tests/pytest.ini']},
     'include_package_data': True
 }


Test patch:
diff --git a/tests/kernels/test_kernels.py b/tests/kernels/test_kernels.py
new file mode 100644
index 000000000..57ed308f6
--- /dev/null
+++ b/tests/kernels/test_kernels.py
@@ -0,0 +1,496 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for the :mod:`pennylane` kernels module.
+"""
+import pennylane as qml
+import pennylane.kernels as kern
+import pytest
+import numpy as np
+from pennylane import numpy as pnp
+import math
+import sys
+
+
+@qml.template
+def _simple_ansatz(x, params):
+    """A simple quantum circuit ansatz."""
+    qml.RX(params[0], wires=[0])
+    qml.RZ(x, wires=[0])
+    qml.RX(params[1], wires=[0])
+
+
+def _mock_kernel(x1, x2, history):
+    """A kernel that memorizes its calls and encodes a fixed values for equal/unequal
+    datapoint pairs."""
+    history.append((x1, x2))
+
+    if x1 == x2:
+        return 1
+    else:
+        return 0.2
+
+
+def _laplace_kernel(x1, x2):
+    """The laplace kernel on scalar data."""
+    return np.exp(-math.fabs(x1 - x2))
+
+
+class TestKernelMatrix:
+    """Tests kernel matrix computations."""
+
+    def test_simple_kernel(self):
+        """Test square_kernel_matrix and kernel_matrix of the _mock_kernel above."""
+        X1 = [0.1, 0.4]
+        X2 = [0.1, 0.3, 0.5]
+
+        K1_expected = pnp.array([[1, 0.2], [0.2, 1]])
+        K2_expected = pnp.array([[1, 0.2, 0.2], [0.2, 0.2, 0.2]])
+
+        K1 = kern.square_kernel_matrix(X1, lambda x1, x2: _mock_kernel(x1, x2, []))
+        K2 = kern.kernel_matrix(X1, X2, lambda x1, x2: _mock_kernel(x1, x2, []))
+
+        assert np.array_equal(K1, K1_expected)
+        assert np.array_equal(K2, K2_expected)
+
+    def test_laplace_kernel(self):
+        """Test square_kernel_matrix and kernel_matrix of the _laplace_kernel above."""
+        X1 = [0.1, 0.4, 0.2]
+        X2 = [0.0, 0.1, 0.3, 0.2]
+
+        K1_expected = pnp.exp(-np.array([[0.0, 0.3, 0.1], [0.3, 0.0, 0.2], [0.1, 0.2, 0.0]]))
+        K2_expected = pnp.exp(
+            -np.array([[0.1, 0.0, 0.2, 0.1], [0.4, 0.3, 0.1, 0.2], [0.2, 0.1, 0.1, 0.0]])
+        )
+
+        K1 = kern.square_kernel_matrix(X1, _laplace_kernel, assume_normalized_kernel=False)
+        K2 = kern.kernel_matrix(X1, X2, _laplace_kernel)
+
+        assert np.allclose(K1, K1_expected)
+        assert np.allclose(K2, K2_expected)
+
+
+class TestKernelPolarity:
+    """Tests kernel methods to compute polarity."""
+
+    def test_correct_calls(self):
+        """Test number and order of calls of the kernel function when computing the
+        polarity, including computation of the diagonal kernel matrix entries."""
+        X = [0.1, 0.4]
+        Y = [1, -1]
+
+        hist = []
+
+        kern.polarity(X, Y, lambda x1, x2: _mock_kernel(x1, x2, hist))
+
+        assert hist == [(0.1, 0.1), (0.1, 0.4), (0.4, 0.4)]
+
+    def test_correct_calls_normalized(self):
+        """Test number and order of calls of the kernel function when computing the
+        polarity, assuming normalized diagonal kernel matrix entries."""
+        X = [0.1, 0.4]
+        Y = [1, -1]
+
+        hist = []
+
+        kern.polarity(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, hist), assume_normalized_kernel=True
+        )
+
+        assert hist == [(0.1, 0.4)]
+
+    def test_polarity_value(self):
+        """Test value of polarity without class label rescaling (1/2)."""
+        X = [0.1, 0.4]
+        Y = [1, -1]
+        pol = kern.polarity(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, []), rescale_class_labels=False
+        )
+        pol_assume = kern.polarity(
+            X,
+            Y,
+            lambda x1, x2: _mock_kernel(x1, x2, []),
+            assume_normalized_kernel=True,
+            rescale_class_labels=False,
+        )
+
+        assert pol == 1.6
+        assert pol == pol_assume
+
+    def test_polarity_value_other_labels(self):
+        """Test value of polarity without class label rescaling (2/2)."""
+        X = [0.1, 0.4]
+        Y = [1, 1]
+        pol = kern.polarity(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, []), rescale_class_labels=False
+        )
+        pol_assume = kern.polarity(
+            X,
+            Y,
+            lambda x1, x2: _mock_kernel(x1, x2, []),
+            assume_normalized_kernel=True,
+            rescale_class_labels=False,
+        )
+
+        assert pol == 2.4
+        assert pol == pol_assume
+
+
+class TestKernelTargetAlignment:
+    """Tests computation of kernel target alignment."""
+
+    def test_correct_calls(self):
+        """Test number and order of calls of the kernel function when computing the
+        kernel target alignment, including computation of the diagonal kernel matrix entries."""
+        X = [0.1, 0.4]
+        Y = [1, -1]
+
+        hist = []
+
+        kern.target_alignment(X, Y, lambda x1, x2: _mock_kernel(x1, x2, hist))
+
+        assert hist == [(0.1, 0.1), (0.1, 0.4), (0.4, 0.4)]
+
+    def test_correct_calls_normalized(self):
+        """Test number and order of calls of the kernel function when computing the
+        kernel target alignment, assuming normalized diagonal kernel matrix entries."""
+        X = [0.1, 0.4]
+        Y = [1, -1]
+
+        hist = []
+
+        kern.target_alignment(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, hist), assume_normalized_kernel=True
+        )
+
+        assert hist == [(0.1, 0.4)]
+
+    def test_alignment_value(self):
+        """Test value of kernel target alignment without class label rescaling (1/3)."""
+        X = [0.1, 0.4]
+        Y = [1, -1]
+
+        alignment = kern.target_alignment(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, []), rescale_class_labels=False
+        )
+        alignment_assume = kern.target_alignment(
+            X,
+            Y,
+            lambda x1, x2: _mock_kernel(x1, x2, []),
+            assume_normalized_kernel=True,
+            rescale_class_labels=False,
+        )
+
+        assert alignment == 1.6 / (2 * math.sqrt(2 + 2 * 0.2 ** 2))
+        assert alignment == alignment_assume
+
+    def test_alignment_value_other_labels(self):
+        """Test value of kernel target alignment without class label rescaling (2/3)."""
+        X = [0.1, 0.4]
+        Y = [1, 1]
+        alignment = kern.target_alignment(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, []), rescale_class_labels=False
+        )
+        alignment_assume = kern.target_alignment(
+            X,
+            Y,
+            lambda x1, x2: _mock_kernel(x1, x2, []),
+            assume_normalized_kernel=True,
+            rescale_class_labels=False,
+        )
+
+        assert alignment == 2.4 / (2 * math.sqrt(2 + 2 * 0.2 ** 2))
+        assert alignment == alignment_assume
+
+    def test_alignment_value_three(self):
+        """Test value of kernel target alignment without class label rescaling
+        on more data (3/3)."""
+        X = [0.1, 0.4, 0.0]
+        Y = [1, -1, 1]
+
+        alignment = kern.target_alignment(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, []), rescale_class_labels=False
+        )
+        alignment_assume = kern.target_alignment(
+            X,
+            Y,
+            lambda x1, x2: _mock_kernel(x1, x2, []),
+            assume_normalized_kernel=True,
+            rescale_class_labels=False,
+        )
+
+        K1 = np.array(
+            [
+                [
+                    1,
+                    0.2,
+                    0.2,
+                ],
+                [0.2, 1, 0.2],
+                [0.2, 0.2, 1],
+            ]
+        )
+        K2 = np.array([[1, -1, 1], [-1, 1, -1], [1, -1, 1]])
+        expected_alignment = np.trace(np.dot(K1, K2)) / math.sqrt(
+            np.trace(np.dot(K1, K1)) * np.trace(np.dot(K2, K2))
+        )
+
+        assert alignment == expected_alignment
+        assert alignment == alignment_assume
+
+    def test_alignment_value_with_normalization(self):
+        """Test value of kernel target alignment with class label rescaling."""
+        X = [0.1, 0.4, 0.0]
+        Y = [1, -1, 1]
+
+        alignment = kern.target_alignment(
+            X, Y, lambda x1, x2: _mock_kernel(x1, x2, []), rescale_class_labels=True
+        )
+        alignment_assume = kern.target_alignment(
+            X,
+            Y,
+            lambda x1, x2: _mock_kernel(x1, x2, []),
+            assume_normalized_kernel=True,
+            rescale_class_labels=True,
+        )
+
+        K1 = np.array(
+            [
+                [
+                    1,
+                    0.2,
+                    0.2,
+                ],
+                [0.2, 1, 0.2],
+                [0.2, 0.2, 1],
+            ]
+        )
+        _Y = np.array([1 / 2, -1, 1 / 2])
+        K2 = np.outer(_Y, _Y)
+        expected_alignment = np.trace(np.dot(K1, K2)) / math.sqrt(
+            np.trace(np.dot(K1, K1)) * np.trace(np.dot(K2, K2))
+        )
+
+        assert alignment == expected_alignment
+        assert alignment == alignment_assume
+
+
+class TestRegularization:
+    """Tests regularization/postprocessing methods."""
+
+    @pytest.mark.parametrize(
+        "input",
+        [
+            (np.diag([1, 0.4])),
+            (np.diag([1, 0.0])),
+            (np.array([[1, -0.5], [-0.5, 1]])),
+        ],
+    )
+    def test_do_nothing_on_non_negative(self, input):
+        """Test thresholding, displacing and flipping matrix to do nothing on PSD matrices."""
+        assert np.allclose(kern.threshold_matrix(input), input)
+        assert np.allclose(kern.displace_matrix(input), input)
+        assert np.allclose(kern.flip_matrix(input), input)
+
+    @pytest.mark.parametrize(
+        "input,expected_output",
+        [
+            (np.diag([1, -1]), np.diag([1, 0])),
+            (np.array([[1, 1], [1, -1]]), np.array([[1.20710678, 0.5], [0.5, 0.20710678]])),
+            (np.array([[0, 1], [1, 0]]), np.array([[1, 1], [1, 1]]) / 2.0),
+        ],
+    )
+    def test_threshold_matrix(self, input, expected_output):
+        """Test thresholding of eigenvalues of a matrix."""
+        assert np.allclose(kern.threshold_matrix(input), expected_output)
+
+    @pytest.mark.parametrize(
+        "input,expected_output",
+        [
+            (np.diag([1, -1]), np.diag([2, 0])),
+            (
+                np.array([[1, 1], [1, -1]]),
+                np.array([[1 + math.sqrt(2), 1], [1, -1 + math.sqrt(2)]]),
+            ),
+            (np.array([[0, 1], [1, 0]]), np.array([[1, 1], [1, 1]])),
+        ],
+    )
+    def test_displace_matrix(self, input, expected_output):
+        """Test displacing/shifting of eigenvalues of a matrix."""
+        assert np.allclose(kern.displace_matrix(input), expected_output)
+
+    @pytest.mark.parametrize(
+        "input,expected_output",
+        [
+            (np.diag([1, -1]), np.diag([1, 1])),
+            (
+                np.array([[1, 1], [1, -1]]),
+                np.array([[math.sqrt(2), 0], [0, math.sqrt(2)]]),
+            ),
+            (np.array([[0, 1], [1, 0]]), np.array([[1, 0], [0, 1]])),
+        ],
+    )
+    def test_flip_matrix(self, input, expected_output):
+        """Test taking the absolute values of eigenvalues of a matrix."""
+        assert np.allclose(kern.flip_matrix(input), expected_output)
+
+    @pytest.mark.parametrize(
+        "input,fix_diagonal,expected_output",
+        [
+            (np.diag([1, -1]), False, np.diag([1, 0])),
+            (
+                np.array([[1, 1], [1, -1]]),
+                False,
+                np.array([[math.sqrt(2) + 1, 1], [1, math.sqrt(2) - 1]]) / 2,
+            ),
+            (np.array([[0, 1], [1, 0]]), False, np.array([[1, 1], [1, 1]]) / 2),
+            (np.diag([1, -1]), True, np.diag([1, 1])),
+            (np.array([[1, 1], [1, -1]]), True, np.array([[1.0, 1.0], [1.0, 1.0]])),
+            # the small perturbation ensures that the solver does not get stuck
+            (np.array([[0, 1.000001], [1, 0]]), True, np.array([[1, 1], [1, 1]])),
+        ],
+    )
+    def test_closest_psd_matrix(self, input, fix_diagonal, expected_output):
+        """Test obtaining the closest positive semi-definite matrix using a semi-definite program."""
+        try:
+            import cvxpy as cp
+
+            output = kern.closest_psd_matrix(input, fix_diagonal=fix_diagonal, feastol=1e-10)
+        except ModuleNotFoundError:
+            pytest.skip(
+                "cvxpy seems to not be installed on the system."
+                "It is required for qml.kernels.closest_psd_matrix"
+                " and can be installed via `pip install cvxpy`."
+            )
+        except cp.error.SolverError:
+            pytest.skip(
+                "The cvxopt solver seems to not be installed on the system."
+                "It is the default solver for qml.kernels.closest_psd_matrix"
+                " and can be installed via `pip install cvxopt`."
+            )
+
+        assert np.allclose(output, expected_output, atol=1e-5)
+
+    @pytest.mark.parametrize(
+        "input",
+        [
+            (np.diag([1, -1])),
+        ],
+    )
+    def test_closest_psd_matrix_import_error(self, input, mocker):
+        """Test import error raising if cvxpy is not installed."""
+        with pytest.raises(ImportError) as import_error:
+            mocker.patch.dict(sys.modules, {"cvxpy": None})
+            output = kern.closest_psd_matrix(input, fix_diagonal=True, feastol=1e-10)
+
+        assert "CVXPY is required" in str(import_error.value)
+
+    @pytest.mark.parametrize(
+        "input,solver",
+        [
+            (np.diag([1, -1]), "I am not a solver"),
+        ],
+    )
+    def test_closest_psd_matrix_solve_error(self, input, solver):
+        """Test verbose error raising if problem.solve crashes."""
+        with pytest.raises(Exception) as solve_error:
+            output = kern.closest_psd_matrix(input, solver=solver, fix_diagonal=True, feastol=1e-10)
+
+        assert "CVXPY solver did not converge." in str(solve_error.value)
+
+
+def depolarize(mat, rates, num_wires, level):
+    """Apply effect of depolarizing noise in circuit to kernel matrix."""
+    if level == "per_circuit":
+        noisy_mat = (1 - rates) * mat + rates / (2 ** num_wires) * np.ones_like(mat)
+    elif level == "per_embedding":
+        noisy_mat = np.copy(mat)
+        for i in range(len(mat)):
+            for j in range(i, len(mat)):
+                rate = rates[i] + rates[j] - rates[i] * rates[j]
+                noisy_mat[i, j] *= 1 - rate
+                noisy_mat[i, j] += rate / (2 ** num_wires)
+
+    return noisy_mat
+
+
+class TestMitigation:
+    """Tests depolarizing noise mitigation techniques."""
+
+    num_wires = 1
+
+    @pytest.mark.parametrize(
+        "input, use_entries, expected_output",
+        [
+            (np.diag([0.9, 0.9]), (0,), np.array([[1, -1 / 8], [-1 / 8, 1]])),
+            (np.diag([0.9, 0.9]), (1,), np.array([[1, -1 / 8], [-1 / 8, 1]])),
+            (np.diag([1.0, 0.9]), None, np.diag([1, 0.9])),
+            (np.diag([1.0, 0.9]), (1,), np.array([[9 / 8, -1 / 8], [-1 / 8, 1.0]])),
+            (
+                depolarize(np.array([[1.0, 0.5], [0.5, 1.0]]), 0.1, num_wires, "per_circuit"),
+                (0,),
+                np.array([[1.0, 0.5], [0.5, 1.0]]),
+            ),
+            (
+                depolarize(np.array([[1.0, 0.5], [0.5, 1.0]]), 0.1, num_wires, "per_circuit"),
+                (1,),
+                np.array([[1.0, 0.5], [0.5, 1.0]]),
+            ),
+        ],
+    )
+    def test_mitigate_depolarizing_noise_single(self, input, use_entries, expected_output):
+        """Test mitigation of depolarizing noise in kernel matrix measuring a single noise rate."""
+        output = kern.mitigate_depolarizing_noise(input, self.num_wires, "single", use_entries)
+        assert np.allclose(output, expected_output)
+
+    @pytest.mark.parametrize(
+        "input, use_entries, expected_output",
+        [
+            (np.diag([0.9, 0.9]), None, np.array([[1, -1 / 8], [-1 / 8, 1]])),
+            (np.diag([1.0, 0.9]), None, np.array([[19 / 18, -1 / 18], [-1 / 18, 17 / 18]])),
+            (np.diag([1.0, 0.9]), (0, 1), np.array([[19 / 18, -1 / 18], [-1 / 18, 17 / 18]])),
+            (
+                depolarize(np.array([[1.0, 0.5], [0.5, 1.0]]), 0.1, num_wires, "per_circuit"),
+                None,
+                np.array([[1.0, 0.5], [0.5, 1.0]]),
+            ),
+        ],
+    )
+    def test_mitigate_depolarizing_noise_average(self, input, use_entries, expected_output):
+        """Test mitigation of depolarizing noise in kernel matrix averaging a single noise rate."""
+        output = kern.mitigate_depolarizing_noise(input, self.num_wires, "average", use_entries)
+        assert np.allclose(output, expected_output)
+
+    @pytest.mark.parametrize(
+        "input, expected_output",
+        [
+            (np.diag([0.9, 0.9]), np.array([[1, -1 / 8], [-1 / 8, 1]])),
+            (np.diag([1.0, 0.9]), np.array([[1, -0.059017], [-0.059017, 1.0]])),
+            (
+                depolarize(np.array([[1.0, 0.5], [0.5, 1.0]]), 0.1, num_wires, "per_circuit"),
+                np.array([[1.0, 0.5], [0.5, 1.0]]),
+            ),
+            (
+                depolarize(
+                    np.array([[1.0, 0.5], [0.5, 1.0]]), [0.2, 0.1], num_wires, "per_embedding"
+                ),
+                np.array([[1.0, 0.5], [0.5, 1.0]]),
+            ),
+        ],
+    )
+    def test_mitigate_depolarizing_noise_split_channel(self, input, expected_output):
+        """Test mitigation of depolarizing noise in kernel matrix estimating individual
+        noise rates per datapoint."""
+        output = kern.mitigate_depolarizing_noise(input, self.num_wires, "split_channel")
+        assert np.allclose(output, expected_output)
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 14344939e..201fa9e35 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -204,7 +204,8 @@ class TestPauliEigs:
     def test_correct_eigenvalues_pauli_kronecker_products_three_qubits(self, pauli_product):
         """Test the paulieigs function for three qubits"""
         assert np.array_equal(
-            pu.pauli_eigs(3), np.diag(np.kron(self.pauliz, np.kron(self.pauliz, self.pauliz)))
+            pu.pauli_eigs(3),
+            np.diag(np.kron(self.pauliz, np.kron(self.pauliz, self.pauliz))),
         )
 
     @pytest.mark.parametrize("depth", list(range(1, 6)))
@@ -239,7 +240,12 @@ class TestArgumentHelpers:
             pass
 
         res = pu._get_default_args(dummy_func)
-        expected = {"c": (2, 8), "d": (3, [0, 0.65]), "e": (4, np.array([4])), "f": (5, None)}
+        expected = {
+            "c": (2, 8),
+            "d": (3, [0, 0.65]),
+            "e": (4, np.array([4])),
+            "f": (5, None),
+        }
 
         assert res == expected
 
@@ -328,7 +334,8 @@ class TestExpand:
     def test_expand_invalid_wires(self):
         """test exception raised if unphysical subsystems provided."""
         with pytest.raises(
-            ValueError, match="Invalid target subsystems provided in 'original_wires' argument"
+            ValueError,
+            match="Invalid target subsystems provided in 'original_wires' argument",
         ):
             pu.expand(U2, [-1, 5], 4)
 
@@ -440,7 +447,8 @@ class TestExpand:
     def test_expand_vector_invalid_wires(self):
         """Test exception raised if unphysical subsystems provided."""
         with pytest.raises(
-            ValueError, match="Invalid target subsystems provided in 'original_wires' argument"
+            ValueError,
+            match="Invalid target subsystems provided in 'original_wires' argument",
         ):
             pu.expand_vector(TestExpand.VECTOR2, [-1, 5], 4)
 
@@ -830,7 +838,8 @@ class TestInv:
     def test_non_operations_in_list(self, arg):
         """Test that the proper error is raised when the argument does not only contain operations."""
         with pytest.raises(
-            ValueError, match="The given operation_list does not only contain Operations"
+            ValueError,
+            match="The given operation_list does not only contain Operations",
         ):
             pu.inv(arg)
 
@@ -842,3 +851,28 @@ class TestInv:
             match=r"Use of qml\.inv\(\) is deprecated and should be replaced with qml\.adjoint\(\)\.",
         ):
             qml.inv(qml.Hadamard(wires=[0]))
+
+
+class TestFrobeniusInnerProduct:
+    @pytest.mark.parametrize(
+        "A,B,normalize,expected",
+        [
+            (np.eye(2), np.eye(2), False, 2.0),
+            (np.eye(2), np.zeros((2, 2)), False, 0.0),
+            (
+                np.array([[1.0, 2.3], [-1.3, 2.4]]),
+                np.array([[0.7, -7.3], [-1.0, -2.9]]),
+                False,
+                -21.75,
+            ),
+            (np.eye(2), np.eye(2), True, 1.0),
+            (
+                np.array([[1.0, 2.3], [-1.3, 2.4]]),
+                np.array([[0.7, -7.3], [-1.0, -2.9]]),
+                True,
+                -0.7381450594,
+            ),
+        ],
+    )
+    def test_frobenius_inner_product(self, A, B, normalize, expected):
+        assert expected == pytest.approx(pu.frobenius_inner_product(A, B, normalize=normalize))