Least squares for a rectangular system, exact-capable. More...

#include <cstddef>
#include <vector>
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"

Include dependency graph for lstsq.h:

Classes
struct	line::LstsqResult< T >
	Outcome of lstsq: the solution and whether the system was rank deficient. More...

Namespaces
namespace	line

Functions
template<class T>
std::vector< std::size_t >	line::rref (Matrix< T > &A, const T &tol)
	Reduced row echelon form of A, in place, returning the pivot columns.
template<class T>
LstsqResult< T >	line::lstsq (const Matrix< T > &A, const std::vector< T > &b, const T &tol)
	Least-squares solution of A x = b, minimum-norm when A is rank deficient.
template<class T>
LstsqResult< T >	line::lstsq (const Matrix< T > &A, const std::vector< T > &b)
	Overload picking the default pivot threshold for the arithmetic in use.

Detailed Description

Least squares for a rectangular system, exact-capable.

MATLAB reaches for qr(A,0) when A has full column rank and for a truncated SVD when it does not. Neither has an exact counterpart: a QR factor contains square roots and a singular value is algebraic, not rational. The two routines here give the SAME vectors by rational means.

Full column rank: the normal equations A'A x = A'b. Their solution IS the least-squares solution, identically, not an approximation of it; the usual objection to them is the squared condition number, which is a floating-point concern and disappears at T = Rational.
Rank deficient: the minimum-norm least-squares solution x = A^+ b, built from a full-rank factorization A = C F obtained by exact elimination,

x = F' (F F')^{-1} (C'C)^{-1} C' b.

For a matrix of rank k this equals the truncated-SVD solution at the same k, so it reproduces MATLAB's fallback rather than approximating it.

RANK DETECTION is where exact and inexact arithmetic genuinely part company. At T = Rational a pivot is either zero or it is not, so the rank is exact and tol is ignored. In floating point a threshold is unavoidable; the default is max(m,n) * 1e-14 * max|A|, in the spirit of MATLAB's max(size(A)) * eps(max(sv)) but computed from the entries rather than from singular values that are not available here.

NOT PORTED, deliberately: MATLAB's callers follow a rank-deficient solve with a RANDOM perturbation retry (rng(23000,'twister'), demands nudged by 1e-10..1e-4 times their scale) and keep whichever perturbed model looks better behaved. That is a floating-point remedy for a floating-point rank test. It has no exact counterpart, it changes the answer, and reproducing it would require MATLAB's Mersenne-Twister double stream bit for bit. The exact pseudoinverse above is what the perturbation is trying to approximate.

Definition in file lstsq.h.

Classes

Namespaces

Functions

Detailed Description