How to express real matrices as linear combinations of unitaries?

Question

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?

Answer 1

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!

Answer 2

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])]
        else:
            # 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
    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()