Number of Givens rotations does not decrease monotonically with increased tolerance

[kevinsung]

Courtesy of @jrm874

import numpy as np
import scipy.linalg

import ffsim


def siam_diag_bath_h1e(n_bath: int) -> np.ndarray:
    hopping = -np.diag(np.ones(n_bath - 1), 1) - np.diag(np.ones(n_bath - 1), -1)
    _, eigvecs = np.linalg.eigh(hopping)
    norb = n_bath + 1
    rot = np.eye(norb)
    rot[1:, 1:] = eigvecs
    h1e = np.zeros((norb, norb))
    h1e[0, 1] = h1e[1, 0] = -1.0  # impurity-bath coupling
    for i in range(1, norb - 1):  # bath chain hopping
        h1e[i, i + 1] = h1e[i + 1, i] = -1.0
    return rot.T @ h1e @ rot


h1e = siam_diag_bath_h1e(n_bath=13)
mat = scipy.linalg.expm(-1j * 0.1 * h1e)

tols = [1e-12, 1e-10, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1]

prev = None
monotonic = True
for tol in tols:
    givens_rotations, _ = ffsim.linalg.givens_decomposition(mat, tol=tol)
    count = len(givens_rotations)
    flag = ""
    if prev is not None and count > prev:
        flag = "  <-- INCREASED (non-monotonic!)"
        monotonic = False
    print(f"tol={tol:<8g} num Givens rotations: {count}{flag}")
    prev = count

print(f"\nMonotonic (non-increasing)? {monotonic}")

tol=1e-12    num Givens rotations: 63
tol=1e-10    num Givens rotations: 73  <-- INCREASED (non-monotonic!)
tol=1e-08    num Givens rotations: 46
tol=1e-07    num Givens rotations: 46
tol=1e-06    num Givens rotations: 55  <-- INCREASED (non-monotonic!)
tol=1e-05    num Givens rotations: 38
tol=0.0001   num Givens rotations: 42  <-- INCREASED (non-monotonic!)
tol=0.001    num Givens rotations: 25
tol=0.01     num Givens rotations: 23
tol=0.1      num Givens rotations: 0

Monotonic (non-increasing)? False