How to express real matrices as linear combinations of unitaries?

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?

• Is $A$ a finite-dimensional matrix? There are $N^2$ linearly independent $N\times N$ unitary matrices, so that should help you. Sep 29 at 14:08
• Yes, A is a finite-dimensional matrix, most probably of the form $2^n * 2^n$ Sep 29 at 14:12
• "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?
– glS
Sep 29 at 14:17
• 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? Sep 29 at 14:25
• 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\}$"
– glS
Sep 29 at 15:41

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!

• 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)
– glS
Sep 29 at 15:44
• 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! Sep 30 at 6:22
• 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. Sep 30 at 10:14
• Thanks! Explained beautifully! Sep 30 at 13:41
• 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. Sep 30 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_width = [(0, 0), (0, matrix.shape[0] - matrix.shape[1])]
else:
pad_width = [(0, matrix.shape[1] - matrix.shape[0]), (0, 0)]
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
else:
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()