Skip to content

nlft_qsp.solvers.convolve_optimize

Convolution optimization algorithm for polynomial completion (taken from github.com/Danimhn/GQSP-Code)

Note

This optional module uses torch (~=2.7.1), which is not included in the package requirements.

Functions:

Name Description
complete

Uses the convolution optimization algorithm to find a complementary polynomial to the given one. See arXiv:2308.01501 for an explanation of the method.

complete(b: Polynomial) -> Polynomial

Uses the convolution optimization algorithm to find a complementary polynomial to the given one. See arXiv:2308.01501 for an explanation of the method.

Parameters:

Name Type Description Default
b Polynomial

The polynomial to complete.

required
Note

Numerical stability is not guaranteed.

Returns:

Type Description
Polynomial

A polynomial \(a(z)\) satisfying \(|a(z)|^2 + |b(z)|^2 = 1\) on the unit circle.

Polynomial

In particular \((a, b)\) will be in the image of the NLFT.

Source code in nlft_qsp/solvers/convolve_optimize.py
def complete(b: Polynomial) -> Polynomial:
    """Uses the convolution optimization algorithm to find a complementary polynomial to the given one. See [arXiv:2308.01501](https://arxiv.org/abs/2308.01501) for an explanation of the method.

    Args:
        b (Polynomial): The polynomial to complete.

    Note:
        Numerical stability is not guaranteed.

    Returns:
        A polynomial $a(z)$ satisfying $|a(z)|^2 + |b(z)|^2 = 1$ on the unit circle.
        In particular $(a, b)$ will be in the image of the NLFT.
    """
    poly = torch.tensor(b.coeffs, dtype=torch.cdouble)

    conv_p_negative = complex_convolve(poly)
    conv_p_negative[poly.shape[0] - 1] -= 1

    # Initializing Q randomly to start with
    initial = torch.randn(poly.shape[0]*2, device=device, requires_grad=True)
    initial = (initial / torch.norm(initial)).clone().detach().requires_grad_(True)

    optimizer = torch.optim.LBFGS([initial], max_iter=1000)

    def closure():
        optimizer.zero_grad()
        loss = objective_torch(initial, conv_p_negative)
        loss.backward()
        return loss

    optimizer.step(closure)

    threshold = closure().item()
    attempts = 0
    while closure().item() > bd.machine_threshold():
        optimizer.step(closure)

        new_thr = closure().item()
        if threshold <= new_thr:
            attempts += 1
            if attempts >= CONV_OPT_MAX_ATTEMPTS:
                break
        else:
            threshold = new_thr
            attempts = 0

    opt_initial = initial.detach()
    real = opt_initial[:len(opt_initial) // 2]
    imag = opt_initial[len(opt_initial) // 2:]

    a = Polynomial(torch.complex(real, imag).tolist())
    a = a.shift(-b.effective_degree())
    a *= np.exp(-1j*np.angle(a[0])) # we want a[0] > 0
    return a