HHL algorithm returns incorrect answer with matrix >= 8*8

[xiang-yu]

Environment

• Qiskit Terra version: 0.19.2 (or 0.20)
• Python version: 3.9.6
• Operating system: MacOS Catalina 10.15.7

What is happening?

When I tried to reproduce the Qiskit HHL tutorial (HHL) with A the 8 by 8 matrix, I got very different HHL solutions compared to the classical one.

Here is the result,
Using the 8 by 8 matrix from the Qiskit HHl tutorial,
HHL solution= [-3.63250174e-04-0.00087696j -7.13722578e-05-0.00017231j
-1.37749833e-04-0.00033256j -3.99987591e-03-0.00965655j
2.21093015e-03+0.00533766j 3.26591171e-03+0.00788461j
-1.77709870e-03-0.0042903j 1.06751932e-03+0.00257722j]
Classical solution= [1.14589783 0.4376935 0.16718266 0.06385449 0.0243808 0.00928793
0.00348297 0.00116099]
Fidelity: 0.007554
Euclidean norm of HHL= 0.015680300674385816 ;Euclidean norm of HHL_XYL= 0.015680300674385823
Euclidean norm of classical solution= 1.2399109034486437

How can we reproduce the issue?

import sys
import os
import numpy as np
from qiskit import QuantumCircuit
from qiskit.algorithms.linear_solvers.hhl import HHL
from qiskit import quantum_info
from qiskit.algorithms.linear_solvers.observables import AbsoluteAverage, MatrixFunctional
from qiskit.quantum_info import state_fidelity
from scipy.sparse import diags

hhl = HHL()

nb = 3
a = 1
b = -1/3
matrix = diags([b, a, b], [-1, 0, 1], shape=(2**nb, 2**nb)).toarray()
vector = np.array([1] + [0]*(2**nb -1))
print('Eigenvalues of the input matix:', np.linalg.eigvalsh(matrix))
num_qubits = int(np.log2(matrix.shape[0]))
HHL_solution = hhl.solve(matrix, vector)
total_qubits = HHL_solution.state.num_qubits
HHL_extract = quantum_info.Statevector(HHL_solution.state).data[2 ** (total_qubits - 1) : 2 ** (total_qubits - 1) + 2 ** num_qubits]
HHL_final = HHL_extract*HHL_solution.euclidean_norm/np.linalg.norm(HHL_extract)
classical_solution = NumPyLinearSolver().solve(matrix, vector / np.linalg.norm(vector))
print('HHL solution', HHL_final)
print('Classical solution', classical_solution.state)
print('Euclidean norm of HHL', HHL_solution.euclidean_norm)
print('Euclidean norm of the classical solution', classical_solution.euclidean_norm')

What should happen?

"HHL_final" and "classical_solution.state" should be very close.
Likewise, the Euclidean from the HHL algorithm and the classical one should be almost the same, i.e., the following commands should return very close results,

print('Euclidean norm of HHL', HHL_solution.euclidean_norm)
print('Euclidean norm of the classical solution', classical_solution.euclidean_norm')

Any suggestions?

@anedumla suggested that ([issue 6880])(#6880) there is a bug for the Qiskit implementation. I quote the suggestion of @anedumla here,
"It seems it has to do with the default Hamiltonian simulation from qiskit. If you use instead

from qiskit.algorithms.linear_solvers import TridiagonalToeplitz
tridi_matrix = TridiagonalToeplitz(nb, a, b)
tridi_solution = HHL().solve(tridi_matrix, vector)

It returns the right solution:

print(tridi_solution.euclidean_norm)

1.227424717239271

The difference between calling HHL with TridiagonalToeplitz or an array is that the former implements an efficient circuit and the latter calls qiskit's unitary circuit, i.e.:

qc = QuantumCircuit(self.num_state_qubits)
evolved = sp.linalg.expm(1j * self.matrix * self.evolution_time)
qc.unitary(evolved, qc.qubits)

So it seems there is a bug somewhere - you could open an issue. Thanks for noticing!"

[anedumla]

I will reply directly here your question from #6880 . In the HHL class the call happens in line 367

 matrix_circuit = NumPyMatrix(matrix, evolution_time=2 * np.pi)

This calls qiskit.algorithms.linear_solvers.matrices.numpy_matrix, and then the call happens in line 216:

qc.unitary(evolved, qc.qubits)

That calls the class qiskit.extensions.unitary. Let me know if you have more questions!

[xiang-yu]

Thanks a lot for pointing me to the problem @anedumla . I have looked at the code but couldn't come up with a solution right away. It would be very helpful if you @anedumla or other Qiskit experts could have a look at this issue.

[xiang-yu]

Could you @anedumla @jakelishman please kindly help have a look at this issue?