I am working on using Variational Quantum Linear Solver (VQLS) for some tasks. Here, we need to represent matrix A as a linear combination of unitaries.

$$ {\bf A} = \Sigma^n_{i=1} c_iA_i $$

My questions are:

  1. Is there any general decomposition method of finding such $A_i$ unitaries? Please note that data is real-valued and not only binary.

  2. Should matrix A be broken based on its basis? (Basis can be eigenbasis of the matrix A or it can be general basis like Pauli matrices)

  3. How do we find the value of n i.e., the number of such unitaries?

  • $\begingroup$ Is $A$ a finite-dimensional matrix? There are $N^2$ linearly independent $N\times N$ unitary matrices, so that should help you. $\endgroup$ Sep 29, 2021 at 14:08
  • $\begingroup$ Yes, A is a finite-dimensional matrix, most probably of the form $ 2^n * 2^n $ $\endgroup$
    – Sajal
    Sep 29, 2021 at 14:12
  • 2
    $\begingroup$ "Should matrix A be broken based on its basis?" I don't understand. Isn't this precisely what are you asking? That is, you are asking how to decompose $A$ with this basis? $\endgroup$
    – glS
    Sep 29, 2021 at 14:17
  • $\begingroup$ No, I am asking to break this matrix in terms of unitary matrices, not necessarily based on basis of the matrix A. I want to know different techniques for performing this action. Also wouldn't breaking based on basis result in vectors rather than matrix? $\endgroup$
    – Sajal
    Sep 29, 2021 at 14:25
  • 1
    $\begingroup$ asking "how to decompose $A$ in terms of a collection of unitaries $\{A_i\}$" to me reads the same as "how to break $A$ based on the unitaries $\{A_i\}$" $\endgroup$
    – glS
    Sep 29, 2021 at 15:41

2 Answers 2


You can select a basis of unitary matrices with respect to which you can decompose your matrix. For example, if your matrix $A$ is $2^n\times 2^n$, then you can select the Pauli basis $$ \sigma_y,\qquad y\in\{0,1,2,3\}^n $$ You can find the decomposition very easily. Notice that if $$ A=\sum_yA_y\sigma_y $$ then calculating $$ \text{Tr}(A\sigma_x)=A_x2^n $$ because $A_x^2=I$ and the traces of all tensor products of Paulis except the all-identities tensor are 0.

Take, as an example, a matrix $$ A=\left(\begin{array}{cccc} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1 \end{array}\right). $$ Here, we have $n=2$, so we take the set of Pauli matrices $$ \Lambda=\{I\otimes I,I\otimes X,I\otimes Y,I\otimes Z,X\otimes I,X\otimes X,X\otimes Y,X\otimes Z,Y\otimes I,Y\otimes Y,Y\otimes Z,Z\otimes I,Z\otimes X,Z\otimes Y,Z\otimes Z\}. $$ For each in turn ($\sigma\in\Lambda$) you just calcuate $\text{Tr}(\sigma A)/4$. For instance, $$ \text{Tr}(X\otimes X\cdot A)=2 $$ so $A_{1,1}=1/2$. Ultimately, you find out that $$ A=\frac{1}{2}(X\otimes X-Y\otimes Y)+\frac14Z\otimes Z+\frac34I\otimes I+\frac{1}{4}(Z\otimes I-I\otimes Z) $$

Of course, if you want to ask a question along the lines of the smallest set of unitaries with which you can decompose a specific $A$, that might be a very different question!

  • $\begingroup$ this paper (draft? writing? pdf?) by Wheeler is probably relevant here: reed.edu/physics/faculty/wheeler/documents/Miscellaneous%20Math/… (sorry, don't know of other non-pdf sources) $\endgroup$
    – glS
    Sep 29, 2021 at 15:44
  • $\begingroup$ I might be asking very trivial question here, sorry for that. But could you shed some light on how you get these $ A_y \sigma_y $. That's really the main point of asking this question. If you can show some example or post a link to where this is happening, it would be a life saver. Thanks! $\endgroup$
    – Sajal
    Sep 30, 2021 at 6:22
  • $\begingroup$ You just find the full set of Pauli operators (there's $4^n$ of them) as the $\sigma_y$. You then find the $A_y$ as I specified - multiply your $A$ by $\sigma_y$ and take the trace. $\endgroup$
    – DaftWullie
    Sep 30, 2021 at 10:14
  • $\begingroup$ Thanks! Explained beautifully! $\endgroup$
    – Sajal
    Sep 30, 2021 at 13:41
  • $\begingroup$ And in case it helps, the only reason the coefficients $A_y$ are real is because the example matrix is Hermitian, not because the matrix elements of $A$ are real. $\endgroup$
    – chrysaor4
    Sep 30, 2021 at 19:26

Adapting the answer above, if anyone needs implementation to this question in Qiskit.

from itertools import product
import numpy as np

from qiskit.quantum_info import Operator
from qiskit.circuit import library

def make_square(matrix):
    # Pad with 0 to make square matrix
    if matrix.shape[0] != matrix.shape[1]:
        if matrix.shape[0] > matrix.shape[1]:
            # Pad columns
            pad_width = [(0, 0), (0, matrix.shape[0] - matrix.shape[1])]
            # Pad rows
            pad_width = [(0, matrix.shape[1] - matrix.shape[0]), (0, 0)]
        matrix = np.pad(matrix, pad_width)    
    return matrix

def get_pauli_bases(dimension):
    pauli = {
        'x': Operator(library.XGate().to_matrix()),
        'y': Operator(library.YGate().to_matrix()),
        'z': Operator(library.ZGate().to_matrix()),
        'i': Operator(library.IGate().to_matrix())

    bases = {}
    # Creates bases matrix in dimension mentioned in parameter
    # For dimension 1, bases = {I, X, Y, Z}
    # For dimension 2, bases = {II, IX, IY, IZ, XX, XY, XZ, YY, YZ, ZZ}
    if dimension == 1:
        return pauli
        for permutation in product(*[list(pauli.keys())]*dimension):
            permutation = "".join(permutation)
            bases[permutation] = pauli[permutation[0]]
            for idx in range(1, len(permutation)):
                bases[permutation] = bases[permutation].tensor(pauli[permutation[idx]])
        return bases

def linear_combination_pauli(matrix):
    matrix = make_square(matrix)
    matrix_len = matrix.shape[0]
    # Assuming matrix dimension is in power of 2
    nqubits = int(np.log2(matrix_len))
    bases = get_pauli_bases(nqubits)
    decomposition = {}
    for base, base_matrix in bases.items():
        decomposition[base] = np.trace(np.dot(base_matrix, matrix)) / matrix_len
    return decomposition, bases

def validate_decomposition(decompositon, bases, original):
    created = sum(coeff*matrix.data for coeff, matrix in zip(decomposition.values(), bases.values()))
    created = np.around(created, 3)
    return (created == original).all()


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.