Can't compile HHL Circuit for 16x16 matrix

[kharazity]

Environment

• Qiskit Terra version: 0.20.1
• Python version: 3.9.7
• Operating system: Ubuntu 21

What is happening?

I am trying to analyze the circuit for HHL for a set of matrix inversion problems. I am currently trying to compile the circuit for solving this 16x16 system, dilated from a 9x9 system. The dilation is performed by padding the off diagonal terms with 0's and putting 1's on the diagonal.

 K = 
[[ 1.9375    +0.j -0.51041667+0.j -0.02083333+0.j -1.02083333+0.j
  -0.51041667+0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [-0.51041667+0.j  0.97916667+0.j -0.51041667+0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [-0.02083333+0.j -0.51041667+0.j  1.9375    +0.j -1.02083333+0.j
   0.        +0.j -0.51041667+0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [-1.02083333+0.j  0.        +0.j -1.02083333+0.j  3.875     +0.j
  -0.02083333+0.j -0.02083333+0.j -1.02083333+0.j -1.02083333+0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [-0.51041667+0.j  0.        +0.j  0.        +0.j -0.02083333+0.j
   0.95833333+0.j  0.        +0.j -0.51041667+0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j -0.51041667+0.j -0.02083333+0.j
   0.        +0.j  0.95833333+0.j  0.        +0.j -0.51041667+0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j -1.02083333+0.j
  -0.51041667+0.j  0.        +0.j  1.9375    +0.j -0.02083333+0.j
  -0.51041667+0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j -1.02083333+0.j
   0.        +0.j -0.51041667+0.j -0.02083333+0.j  1.9375    +0.j
  -0.51041667+0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j -0.51041667+0.j -0.51041667+0.j
   0.97916667+0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  1.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  1.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  1.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   1.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  1.        +0.j  0.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  1.        +0.j  0.        +0.j]
 [ 0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j
   0.        +0.j  0.        +0.j  0.        +0.j  1.        +0.j]]

with forcing term:

f_vec =  [0.87758256 0.54030231 0.47415988 0.77015115 1.         0.29192658
 0.87758256 0.47415988 0.54030231 0.         0.         0.
 0.         0.         0.         0.        ]

How can we reproduce the issue?

You can reproduce the issue by

import qiskit as qk
from qiskit.algorithms.linear_solvers.hhl import HHL
tol = .1
q_circ = HHL(epsilon=tol).construct_circuit(K, f_vec)

with K and f_vec defined as above. It takes a really long time to just compile the program. This seems to be unreasonable, it shouldn't intend too great an effort to spit out a qasm file, right?

What should happen?

I would expect that compiling the circuit shouldn't be a major overhead. Simulation, sure, but compilation ought not to be.

Any suggestions?

You can compile HHL sequentially. The QPE routine is essentially repeated twice, so you should only have to compile that part of the circuit once and then apply the Hermitian conjugate QPE after the controlled rotation step. I'm not sure where the major bottleneck is.

[woodsp-ibm]

@anedumla Any comment/thought?

[kharazity]

I also wanted to add some relevant details, the condition number of the above matrix is ~46 and the norm of f_vec ~2. Epsilon, 1e-1. The number of qubits is certainly >4, because of the encoding qubits used to store the values for the controlled rotations. I've been using the log(kappa^2/eps) as a rule of thumb for estimating the number of encoding qubits, so in this case ~4 encoding qubits, 1 ancilla, and 4 qubits for storing the 16 states of the matrix. In all, I'd anticipate about 10 qubits are involved in this computation. Does that volume cause any difficulty?

[kharazity]

Just bumping this thread to see if anyone has any suggestions here. A complete circuit description isn't even that important, if one could just even cost out the number of 2-qubit unitaries that would be useful.

[anedumla]

The problem seems to be QPE with the unitary returned by QuantumCircuit.unitary(exp(k)), if instead of your K you create some 16x16 TridiagonalToeplitz matrix the circuit gets constructed very fast.

[kharazity]

I see. That won't really work for what I'm doing, but thank you for letting me know!

[kharazity]

I started just playing with the ToeplitzTridiagonal for HHL as well. I wanted to compile it to qasm, but get the following error:

~/.local/lib/python3.9/site-packages/qiskit/circuit/quantumcircuit.py in qasm(self, formatted, filename, encoding)
   1675 
   1676         # insert gate definitions
-> 1677         string_temp = _insert_composite_gate_definition_qasm(
   1678             string_temp, existing_composite_circuits, self.extension_lib
   1679         )

~/.local/lib/python3.9/site-packages/qiskit/circuit/quantumcircuit.py in _insert_composite_gate_definition_qasm(string_temp, existing_composite_circuits, extension_lib)
   4781             qasm_string = instruction._qasm_definition
   4782         else:
-> 4783             qasm_string = _get_composite_circuit_qasm_from_instruction(instruction)
   4784         gate_definition_string += "\n" + qasm_string
   4785 

~/.local/lib/python3.9/site-packages/qiskit/circuit/quantumcircuit.py in _get_composite_circuit_qasm_from_instruction(instruction)
   4747             ["q%i" % index for index in [definition_bit_labels[qubit] for qubit in qargs]]
   4748         )
-> 4749         composite_circuit_gates += f"{sub_instruction.qasm()} {gate_qargs}; "
   4750 
   4751     if composite_circuit_gates:

~/.local/lib/python3.9/site-packages/qiskit/circuit/instruction.py in qasm(self)
    451             name_param = "{}({})".format(
    452                 name_param,
--> 453                 ",".join([pi_check(i, ndigits=8, output="qasm") for i in self.params]),
    454             )
    455 

~/.local/lib/python3.9/site-packages/qiskit/circuit/instruction.py in <listcomp>(.0)
    451             name_param = "{}({})".format(
    452                 name_param,
--> 453                 ",".join([pi_check(i, ndigits=8, output="qasm") for i in self.params]),
    454             )
    455 

~/.local/lib/python3.9/site-packages/qiskit/circuit/tools/pi_check.py in pi_check(inpt, eps, output, ndigits)
    164         return str_out
    165 
--> 166     complex_inpt = complex(inpt)
    167     real, imag = map(normalize, [complex_inpt.real, complex_inpt.imag])
    168 

TypeError: only length-1 arrays can be converted to Python scalars

Here's the code:

from qiskit.algorithms.linear_solvers import HHL
from qiskit.algorithms.linear_solvers import TridiagonalToeplitz
import numpy as np

Nq = 4
mat = TridiagonalToeplitz(Nq, -2, 1)
b_vec = np.pi*np.ones(2**Nq)
circ = HHL().construct_circuit(mat,b_vec)
circ.qasm()

[rht]

I did a line profiling on HHL.construct_circuit on 3 qubits. The most time-consuming parts are indeed the PhaseEstimation construction, and then its inverse:

   459                                                   # QPE
   460         1    7797279.0 7797279.0     49.8          phase_estimation = PhaseEstimation(nl, matrix_circuit)
   461         1          8.0      8.0      0.0          if na > 0:
   462                                                       qc.append(phase_estimation, ql[:] + qb[:] + qa[: matrix_circuit.num_ancillas])
   463                                                   else:
   464         1     544355.0 544355.0      3.5              qc.append(phase_estimation, ql[:] + qb[:])
   465                                                   # Conditioned rotation
   466         1         17.0     17.0      0.0          if self._exact_reciprocal:
   467         1     165837.0 165837.0      1.1              qc.append(reciprocal_circuit, ql[::-1] + [qf[0]])
   468                                                   else:
   469                                                       qc.append(
   470                                                           reciprocal_circuit.to_instruction(),
   471                                                           ql[:] + [qf[0]] + qa[: reciprocal_circuit.num_ancillas],
   472                                                       )
   473                                                   # QPE inverse
   474         1          8.0      8.0      0.0          if na > 0:
   475                                                       qc.append(phase_estimation.inverse(), ql[:] + qb[:] + qa[: matrix_circuit.num_ancillas])
   476                                                   else:
   477         1    6621489.0 6621489.0     42.3              qc.append(phase_estimation.inverse(), ql[:] + qb[:])

However, the most time-consuming part of the PhaseEstimation construction is actually the .control() step.
Here is the line profile of https://github.com/Qiskit/qiskit-terra/blob/e4e4646f18d6e7d5dae0bc5308bcb5f5c9aa206d/qiskit/circuit/library/phase_estimation.py#L94

    94         6         36.0      6.0      0.0          for j in range(num_evaluation_qubits):  # controlled powers
    95         5      33022.0   6604.4      0.2              a = unitary.power(2**j)
    96         5    7257376.0 1451475.2     45.7              a.control()
    97         5    7973177.0 1594635.4     50.2              circuit.compose(unitary.power(2**j).control(), qubits=[j] + qr_state[:], inplace=True)

[rht]

I think I can get 280x speedup of that last line by replacing

circuit.compose(unitary.power(2**j).control(), qubits=[j] + qr_state[:], inplace=True)

with

evolved = sp.linalg.expm(1j * unitary.matrix * unitary.evolution_time)
out = evolved.copy()
for i in range(2**j - 1):
    out *= evolved
out = control_unitary(out)
circuit.unitary(out, [j] + qr_state[:])

where control_unitary is defined as

        def control_unitary(u):
            identity = np.eye(len(u))
            zerozero = [[1.0, 0.0], [0.0, 0.0]]
            oneone = [[0.0, 0.0], [0.0, 1.0]]
            return np.kron(identity, zerozero) + np.kron(u, oneone)

I tested the correctness of the result. For 2 qubits, the first 2 elements of my statevector matches, but for 3 qubits, the elements of my statevector are so small (~1e-16) that even though they are different, it doesn't make sense because they are effectively 0.0. And so, not sure.