HHL - the result is correct for one matrix but wrong for another one

Question

I tried to run HHL algorithm in new version of Qiskit (0.29). Firstly I tried to solve a diagonal system 4x4 with matrix [[1,0,0,0],[0,2,0,0],[0,0,3,0],[0,0,0,4]] and right side [1,4,0.5,3]. Clearly, the solution is [1, 2, 0.1667, 0.75]. In this case, everything was according to the expectation. Then I switched to matrix [[-1,2,0,0],[2,-1,2,0],[0,2,-1,2],[0,0,2,-1]] with right side [1,1,1,1]. Expected solution is [0.2, 0.6, 0.6, 0.2], however [0.0988, 0.1083, 0.119, 0.0884] is returned.

There is no error during a code execution.

Here is my code:

from qiskit import ClassicalRegister, QuantumRegister, QuantumCircuit, Aer, execute, IBMQ
from qiskit.compiler import transpile, assemble
from qiskit.quantum_info import Statevector

from qiskit.algorithms.linear_solvers.hhl import HHL

import numpy as np
import math as m

def dec2bin(decimal, bits): #auxiliary function for conversion dec to bin
    b = bin(decimal)
    b = b[2:len(b)]
    while len(b) < bits:
        b = '0' + b
    return b

processor = Aer.backends(name='qasm_simulator')[0]

#1st linear system
matrix_A = np.array([[1,0,0,0],[0,2,0,0],[0,0,3,0],[0,0,0,4]])
vector_b = np.array([1,4,0.5,3]) #x = [1, 2, 0.1667, 0.75]

#2nd linear system
#matrix_A = np.array([[-1,2,0,0],[2,-1,2,0],[0,2,-1,2],[0,0,2,-1]])
#vector_b = np.array([1,1,1,1])  #x = [0.2, 0.6, 0.6, 0.2]


n = len(matrix_A[2,:]) #number of variables (assuming n=2^k)
qubitsForResults = m.ceil(m.log(n,2)) #number of necessary qubits

bNorm = np.linalg.norm(vector_b) #norm of right-side

hhlSolver = HHL(quantum_instance = processor)
hhlCircuit = hhlSolver.construct_circuit(matrix_A, vector_b) #preparing circuit in general form

circuitToRun = transpile(hhlCircuit,basis_gates=['id', 'rz', 'sx', 'x', 'cx']) #adapt circuit to real quantum processor
circuitToRun.measure_all() #add measurement

#execute the circuit
Shots = 8192
M = execute(circuitToRun, processor, shots = Shots).result().get_counts(circuitToRun)

usedQubits = circuitToRun.width()/2 #width - num of qubits and classical bit summed
zeros = dec2bin(0,usedQubits - qubitsForResults - 1) #working qubits

for i in range(0,n):
    #ancilla qubit has to be |1>, working qubits |0>, results are in last "qubitsForResults" qubits
    #e.g. if qubitsForResults and total qubits is 6, then 3 are working qubits
    #results are stored in |1 000 00>, |1 000 01>, |1 000 10>, |1 000 11>
    indx = '1' + zeros + dec2bin(i,qubitsForResults)
    #bNorm*m.sqrt(M[indx]/Shots) - we need probability amplitudes => sqrt
    #and adjust them by norm of vector b to get "unnormalized" solution x
    print('x_' + str(i) + ': ' + str(round(bNorm*m.sqrt(M[indx]/Shots),4)))

Note that my code works only in case the solution of the system is real and non-negative.

Could anybody help me to understand where is a mistake?

---

EDIT: I also tried to use matrix [[1,0,0,0],[0,2,0,0],[0,0,3.5,-0.5],[0,0,-0.5,3.5]] with right side [1 1 1 1], expected solution is [1 0.5 0.333 0.333]. In this case the code works fine. I prepared this matrix by a unitary transformation of $\text{diag}(1,2,3,4)$. As a result, eigenvalues of the matrix are 1, 2, 3 and 4. So, it seems that the code works correctly if the matrix of the linear system has interger eigenvalues. I observed similar behavior also for 2x2 matrices.

Is it possible that the errors are caused by lack of qubits used for eigenvalue approximation? In case of integer eigenvalues, there is no loss of accuracy if we have enough qubits. In this case the max. eigenvalue is 4 (100b), so three working qubits are sufficient.

---

EDIT 2: I experimented further with 2x2 matrices and directly with state vector instead of probabilities. I found out that (note that right side is [1, 1.2]):

• if matrix is [[1.5, 0.5],[0.5, 1.5]] then HHL returns correct result, note that eigenvalues of the matrix are 1 and 2
• for matrix [[2.5, 1.5],[1.5, 2.5]] everything is also fine, eigenvalues are 4 and 1
• for matrix [[2/3, -1/3],[-1/3, 2/3]] the results are wrong, eigenvalues are 1/3 and 1
• finally, for matrix [[-0.5, -1.5],[-1.5, -0.5]] the result is (0.45, 0.65) but expected is (-0.65, -0.45), eigenvalues of the matrix are -2 and 1. Result (0.45, 0.65) is same as in case of matrix from the first bullet point with eigenvalues 1 and 2

To conclude, it seems that HHL works well only for matrices with positive integer eigenvalues.

I also had a look at the circuits prepared with HHL method and realized that number of working qubits is equivalent to $\log_2(\max\{|\lambda_i|\})$. There is no qubit for sign or additional qubits to represent fractional eigenvalues.

So, am I right that current HHL implementation in Qiskit is intended only for matrices with positive integer eigenvalues?

Answer 1

Based on comment by Lena, I found out that HHL implementation in Qiskit cannot be used for matrices with negative eigenvalues for time being. As stated here

The current version of HHL doesn't support matrices with negative eigenvalues, which seems to be the reason this example is failing. That's a feature we're looking to integrate and will soon be working on!

During experiments, I also realized that current Qiskit HHL implementation has sometimes difficulties with matrices having non-integer eigenvalues. So, this seems to a matter of accuracy and low number of available qubits.

To conclude, my code works well, however, I used it in cases the current HHL implementation is not suitable for.