HHL algorithm returns incorrect answer with matrix >= 8*8 #7926
Description
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!"