num::inverse_norm1_estimate

Defined in header <linear/condition.hpp>

template <class F, std::floating_point T = real>
requires requires(const F &factor, const vec<T> &b, vec<T> &x) {
  solve(factor, b, x);
  solve_transpose(factor, b, x);
}
T inverse_norm1_estimate(const F &factor, idx n);

Estimate \(\|A^{-1}\|_1\) from a factorization of the n-by-n matrix \(A\).

Hager's method with Higham's refinements, as LAPACK's dlacn2: \(\|A^{-1}\|_1\) is the maximum of \(\|A^{-1}x\|_1\) over the unit 1-norm ball, a convex function maximized at a vertex \(e_j\). Each step solves with \(A\), then with \(A^T\) to find the coordinate \(j\) along which that norm grows fastest, and stops when no vertex improves. That costs a handful of solves, \(O(n^2)\) for a dense factor, instead of the \(O(n^3)\) of forming \(A^{-1}\). The result is a lower bound, usually exact or within a factor of 3.

T is the precision of the solves, which lets a float factorization be measured too.