VQLS Algorithm for linear solver

NicoRenaud commented 2 years ago

What are the advantages/drawback compared to HHL ? Any use in the calibration pipeline ?

References

https://qiskit.org/textbook/ch-paper-implementations/vqls.html
Variational Quantum Linear Solver. Arxiv
Decent explanation of the Variational Principle in Qiskit.
- An even better (and shorter) one is available on Sakurai's Modern Quantum Mechanics, episode 5.4.

PabRod commented 2 years ago

The purpose is to solve the linear system for x, for a given A and b:

$A \left|{x}\right\rangle = \left|{b}\right\rangle$

The key idea is to find the ground energy for the projection hamiltonian:

$H = 1 - \left|{b}\right\rangle \left\langle{b}\right|$

Instead of trying to solve directly, we'd rather apply the variational principle to the functional:

$\epsilon[\phi] \equiv \frac{\left\langle{\phi}\right|H\left|{\phi}\right\rangle}{\left\langle{\phi}\right|\left|{\phi}\right\rangle}$

Where:

$\left|{\phi}\right\rangle \equiv A \left|{\psi(k)}\right\rangle$ $\left|{\psi(k)}\right\rangle \equiv V(k) \left|{0}\right\rangle$

and

$\left|{b}\right\rangle \equiv U \left|{0}\right\rangle$

Additionally, it seems to be important to express A as a linear combination of unitary matrices. Nothing is mentioned about the unicity nor existence of this decomposition.

NicoRenaud commented 2 years ago

Msc thesis that applies VQLS to FEM and spend some time on matrix decomposition :

Publicly available here.

PabRod commented 2 years ago

Decomposition in unitary matrices

According to this, the decomposition in unitaries can always be done using 4 matrices or less. They also explain the algorithm.

Quick and dirty implementation

def unitary_decomposition(X):
    """ Unitary decomposition

    Decomposes the complex normalized matrix X into four unitary matrices
    """
    def get_real(x):
        """ Get the real part of x"""
        return 0.5 * (x + np.conjugate(x))

    def get_imag(x):
        """ Get the imaginary part of x """
        return 0.5/(1j) * (x - np.conjugate(x))

    def aux(x):
        """ Auxiliary function

        Performs a matrix operation that we'll need later
        """
        I = np.eye(len(x))
        return 1j*scipy.linalg.sqrtm(I - x**2)

    # Normalize
    norm = np.linalg.norm(X)
    X_n = X / norm

    # Apply the algorithm as described in
    # https://math.stackexchange.com/questions/1710247/every-matrix-can-be-written-as-a-sum-of-unitary-matrices/1710390#1710390
    B = get_real(X_n)
    C = get_imag(X_n)

    ## Get the matrices
    UB = B + aux(B)
    UC = C + aux(C)
    VB = B - aux(B)
    VC = C - aux(C)

    ## Get the coefficients
    cb = norm * 0.5
    cc = cb * 1j

    ## Return
    return [[[cb], UB], 
            [[cb], VB],
            [[cc], UC],
            [[cc], VC]]

def unitary_recomposition(decomposed):
    """ Rebuilds the original matrix from the decomposed one """
    reference = decomposed[0][1]

    temp = np.zeros(len(reference))
    for line in result:
        coef = line[0]
        mat = line[1]
        temp = temp + coef * mat

    return temp

NicoRenaud commented 2 years ago

It's done and it's part of qalcore. We might want to upstream it to qiskit

NLESC-quantum / quantum_comp