🐛 Numerical instability in DD package for large-scale systems

[burgholzer]

Environment information

Dates back to the very origin of the (QM)DD package.

Description

The code down below demonstrates an instability in the matrix DD normalization.

The script just creates a bunch of Hadamards and then reverses them.
Up until (including) 81 qubits, this produces results as expected. The result is equal to the identity and the top edge weight is one.
For n=81+k, the top edge weight will be sqrt(2)^k ... Meaning that, for example, at n=128 the top edge weight is equal to roughly a million 😅
The main reason for this is that the first bunch of H's generates a DD that has a top edge weight of (1/sqrt(2))^n until n hits 81 where the default numerical tolerance kicks in and the value remains constant at 6.4311e-13, which is just above the default tolerance. Now each qubit above 81 in the second series of H's induces an error of sqrt(2).
In extreme cases, this may even overflow the range of double. Particularly, for n>=2129, sqrt(2)^(n-81) overflows the largest possible double precision floating point value (1.7976931348623158e+308).

Possible consequences of the above include

• NaN values after overflows
• DDs collapsing to the all-zero DD (because the top-weight underflows the tolerance) and subsequent segfaults
• Performance problems due to polluting the 1.0 unique table bucket.

Expected behavior

Matrix normalization should be stable and not induce these kinds of errors.

This will most surely require the development of a new normalization scheme for matrix DDs (probably inspired by the one currently used for vectors), which, from past experience, might not be that easy to get right.

How to Reproduce

Run the following test

TEST(DDPackageTest, MatrixNormalizationRegression) {
  const qc::Qubit nqubits = 100U;
  for (auto q = 1U; q <= nqubits; ++q) {
    auto dd = std::make_unique<dd::Package<>>(q);
    auto f = dd->makeIdent();
    for (qc::Qubit i = 0; i < q; ++i) {
      const auto h = dd->makeGateDD(dd::H_MAT, i);
      f = dd->multiply(f, h);
    }
    for (qc::Qubit i = 0; i < q; ++i) {
      const auto h = dd->makeGateDD(dd::H_MAT, q - i - 1);
      f = dd->multiply(f, h);
    }
    std::cout << "Qubits: " << q << ", ";
    std::cout << "Top edge weight: " << f.w.toString(false) << ", ";
    std::cout << "Is identity: " << f.isIdentity(false) << "\n";
    const auto w = static_cast<dd::ComplexValue>(f.w);
    EXPECT_TRUE(w.approximatelyEquals({1., 0.}));
    EXPECT_TRUE(f.isIdentity(false));
  }
}

[burgholzer]

It turns out that this issue is not limited to the matrix DD case and its not solely the matrix normalization that is at fault.
A similar phenomenon also happens when conducting simulations of an even simpler circuit, i.e.

OPENQASM 2.0;
include "qelib1.inc";

qreg q[81];

h q;

works as expected, but, starting from $81+k$ qubits with $k>0$, the top edge weight starts to grow exponentially with a value of $\sqrt{2}^k$.

The main cause for this issue is the tolerance that is being employed in the package.
The tolerance defaults to $2^{-52} * 1024 = 2^{-42} = 2.274e^{-13}$.
Now,

• $\frac{1}{\sqrt{2}}^{81} = 6.431e^{-13}$,
• $\frac{1}{\sqrt{2}}^{82} = 4.547e^{-13} = 2^{-41}$,
• $\frac{1}{\sqrt{2}}^{83} = 3.216e^{-13}$,
• $\frac{1}{\sqrt{2}}^{84} = 2.274e^{-13} = 2^{-42}$.

So, at $84$ qubits, each individual amplitude after applying the Hadamard transform is as small as the tolerance itself. Such a number would be considered approximatelyZero, which probably leads to all kinds of errors. These errors are likely to be more pronounced in the matrix-only case, given that the normalization scheme explicitly pulls out the common factor of $\frac{1}{\sqrt{2}}$ from each Hadamard gate and makes that the new top edge-weight.

Two questions that pop up:

• Why does this also affect the vector decision diagrams where the normalization scheme does not propagate the common factor through the decision diagram?
• Why does the error already show at $82$ qubits and not only at $84$ qubits and beyond?

I'll try to address them in a follow-up comment.

[burgholzer]

Ok. Just noticed, that this might, in fact, be only due to the matrix normalization.
The following circuit works smoothly

OPENQASM 2.0;
include "qelib1.inc";

qreg q[81];
qreg k[1];

h q;
h k;

The reason for the circuit from the previous comment triggering the issue is that the statement

h q;

Is translated to a CompoundOperation in the resulting QuantumComputation.
During simulation getDD(compOp,...) builds the full matrix DD for the compound operation before applying it to the state.
That means, if an error is introduced in the MxM multiplications that form the overall DD for the Hadamard transform, this error is propagated to the state as a result of the subsequent MxV multiplication.

That eliminates the first question and reduces the scope for the error hunt.

[burgholzer]

As for why the error already shows at $82$ qubits:
$$\frac{1}{\sqrt{2}}^{81} - \frac{1}{\sqrt{2}}^{82} = 1.884e^{-13} < 2.274e^{-13} = 2^{-42} = \varepsilon$$
Thus, the lookup of the top edge weight for the matrix DD after applying the $82^{nd}$ Hadamard gate onward will just return $\frac{1}{\sqrt{2}}^{81}$ and, hence, incur an error of $\frac{1}{\sqrt{2}}$.

This error could, in principle, be fixed by adopting a new normalization scheme that is similar to the vector normalization scheme and does not propagate the greatest common divisor upward.

Even then, one has to be very careful as, during DD addition, all edge weights along a path to a terminal are multiplied together before the addition is carried out and the resulting edge weight is put through normalization and a unique table lookup.
Extra care has to be taken that the multiplication of the edge weights as well as the addition itself is performed with higher numerical accuracy without interfering via the tolerance.
Anyone investigating this issue should carefully consider the CachedEdge normalization routines and make sure that these operations happen with the highest possible accuracy.