bertini.linalg.solve / PartialPivLU: non-square mp input corrupts the heap instead of raising

[ofloveandhate]

Summary

bertini.linalg.solve (and the exposed PartialPivLU) accept non-square multiprecision matrices and hand them straight to Eigen::PartialPivLU, which requires a square invertible matrix. Depending on the shape this either corrupts the heap (free(): invalid pointer, SIGABRT — seen as a bare SIGSEGV in a larger program) or silently returns a wrong answer.

The double path is unaffected: it routes through numpy and raises LinAlgError correctly. So the mp and double paths of the same public API disagree about what is an error.

A crash is the worst outcome here because it is uncatchable from Python — a try: ... except Exception: around the call cannot save you, so a single bad shape kills the whole process with no traceback.

Minimal reproducer

import numpy as np, bertini, bertini.linalg as linalg
from bertini.multiprec import complex_mp as C

bertini.default_precision(30)

A = np.array([[C('2e-6'), C('-3e-6'), C('-1e-12')]])   # 1 x 3
b = np.array([C('5e-13')])                             # length 1

linalg.solve(A, b)      # -> free(): invalid pointer  (SIGABRT)

Control: the same call with double raises properly.

linalg.solve(np.array([[2e-6, -3e-6, -1e-12]]), np.array([5e-13]))
# LinAlgError: Last 2 dimensions of the array must be square

Observed behaviour by shape

complex_mp, well-conditioned entries, each case run in its own subprocess:

shape	linalg.solve	linalg.lu(A).solve(b)	double (either)
3 x 3	ok	ok	ok
1 x 3	CRASH	silently wrong answer	raises LinAlgError
3 x 1	silently wrong answer	CRASH	raises LinAlgError
2 x 4	CRASH	silently wrong answer	raises LinAlgError

lstsq, qr and svd handle all of these shapes correctly (they are backed by
ColPivHouseholderQR / JacobiSVD, which are defined for rectangular input).

Cause

python_bindings/src/linalg_export.cpp:

template<typename T>
Vec<T> SolveLinearSystem(Mat<T> const& A, Vec<T> const& b)
{
    return Vec<T>(A.partialPivLu().solve(b));   // no shape check
}

Eigen::PartialPivLU is documented as requiring a square invertible matrix and
eigen_asserts it — but assertions are compiled out under NDEBUG in a release build,
so a non-square argument is undefined behaviour rather than a caught precondition. The
same applies to eigenpy::PartialPivLUSolverVisitor<...>::expose("PartialPivLU") a few
lines below, which is what the lu() path reaches.

Why it bites in practice

Non-square systems legitimately flow around b2 — MakeMovingHomotopy accepts them by
design (#258) — so "just don't do that" is not much of a guard. The way I hit it: a
damped-Newton loop that calls linalg.solve(J, residual) inside try/except Exception
was handed the 1x3 Jacobian of a 1-function/3-variable system. The except cannot catch
heap corruption, so the process died with no traceback and nothing in the log.

Suggested fix

1. Shape-check in SolveLinearSystem before touching Eigen — raise a Python exception
   (std::invalid_argument / ValueError) naming the shape, e.g. "solve expects a
   square matrix; got 1 x 3 — use lstsq for rectangular systems".
2. Same guard for the exposed PartialPivLU / PartialPivLUReal visitors (a thin
   wrapper around the constructor, or a documented precondition plus a check).
3. Worth auditing the other exposed decompositions for preconditions that are only
   eigen_asserted, since those all vanish under NDEBUG.

More generally: it should not be possible to segfault the library from Python at all.

Environment

• b2 VERSION 3.5.0.dev0, branch feature/eval-precision-tolerant @ c24bdc2e
• installed dist bertini2 2.0.2
• eigenpy 3.12.0
• Python 3.14.5, Linux aarch64