Typed dense decompositions and solvers

Status: stable in 1.6.0. All implementations are native Object Pascal and need no external numerical library. Storage, views, kernels, and conversion rules remain documented in Typed dense matrices.

60-second least-squares solve

uses

  AlgebraLib.DenseMatrices,

  AlgebraLib.DenseDecompositions;



var

  A, B, X: IDenseDoubleMatrix;

  Info: TDenseSolveDiagnostics;

begin

  A := TDenseDoubleMatrix.FromValues(4, 2,

    [1.0, 1.0, 1.0, 2.0, 1.0, 3.0, 1.0, 4.0]);

  B := TDenseDoubleMatrix.FromValues(4, 1, [6.0, 5.0, 7.0, 10.0]);

  X := LeastSquares(A, B, Info);

  WriteLn('intercept=', X[0, 0]:0:3, ' slope=', X[1, 0]:0:3);

  WriteLn('residual=', Info.ResidualNorm:0:6);

end.

Expected output:

intercept=3.500 slope=1.400

residual=2.049390

This fits a line without forming normal equations or an inverse. The complete solver-selection example also handles a rank-deficient calibration with SVD and interprets a covariance eigensystem.

Choose a dense solver

Do not choose SVD merely because it is the most general option. QR is the predictable default for a full-rank tall least-squares problem. Choose CPQR when rank must be diagnosed, and SVD when the solution must be minimum norm or rank deficiency is part of the problem.

Sparse, iterative, and matrix-free paths are stable in 2.2.0, as are partial eigensystems and the Hessenberg, Schur, real nonsymmetric, and generalized spectral APIs described below. LDLT, update/downdate, public workspace, parallel/SIMD, GPU, and external BLAS/LAPACK paths remain unsupported. The compatibility IMatrix methods do not change the typed API contracts.

Common contracts

All indices and shapes use zero-based SizeInt. Inputs must be finite. Factorisation clones its source; later source mutation cannot alter the factor. Accessors return copies, and solves never overwrite coefficients or right-hand sides. A factor has no mutable last-result state, so concurrent reads and solves are safe when callers do not concurrently mutate a shared input matrix.

Convenience calls allocate and factor once. A reusable factor retains O(mn) or O(n²) private storage and avoids repeated O(n³) work, but each solve returns a newly allocated result. The 1.6 API has no destructive factorisation or public workspace form.

Every decomposition supports matching Single, Double, TSingleComplex, and TComplex paths wherever the mathematics applies. A typed factor never converts through nested Double storage or IMatrix. Real symmetric eigensystems apply to the real paths; Hermitian eigensystems apply to the complex paths.

Invalid nil/shape/option/tolerance input, non-finite entries, singular factors, invalid structure, and convergence failure raise EDenseMatrixError before a result is returned. User tolerances must be finite and nonnegative; the default is selected by a negative value and scales with precision, dimensions, and the leading pivot/diagonal/singular value.

Diagnostics

TDenseSolveDiagnostics contains:

Method: string;

NumericalRank: SizeInt;

IsRankDeficient: Boolean;

Tolerance: Double;

ConditionIndicator: Double;

ResidualNorm: Double;

BackwardError: Double;

ResidualNorm is ||B-A*X||F. BackwardError is the scale-aware normalized quantity

||B-A*X||F / (||A||F*||X||F + ||B||F).

ConditionIndicator is the smallest accepted to largest relevant pivot, factor diagonal, or singular value. A smaller value warns of greater sensitivity, but it is deliberately not presented as a formal reciprocal condition-number estimate. No diagnostic forms an inverse.

LU, Cholesky, QR/CPQR, and SVD reusable factors expose SolveWithInfo, SolveLeastSquaresWithInfo, or SolveMinimumNormWithInfo. The matching methods without WithInfo return the same solution while omitting residual calculation.

Triangular solve reference

TDenseTriangle is dtLower or dtUpper; TDenseDiagonal is ddNonUnit or ddUnit; TDenseTranspose is dtNoTranspose, dtTranspose, or dtConjugateTranspose.

SolveTriangular(A,B,Triangle,Transpose,Diagonal) accepts square n x n A and n x r B, including n=0 or r=0. Entries outside the selected triangle are ignored. Real transpose and conjugate transpose have the same mathematics. Time is O(n²r); the result needs O(nr) additional storage.

LU and Cholesky reference

The four LU handles are IDenseSingleLU, IDenseDoubleLU, IDenseSingleComplexLU, and IDenseComplexLU. They expose Size, copied L, U, and Permutation, PivotRatio, ConditionIndicator, IsIllConditioned, Solve, and SolveWithInfo. The permutation convention is P*A=L*U; Permutation[i] is the source row occupying factor row i.

The four Cholesky handles are IDenseSingleCholesky, IDenseDoubleCholesky, IDenseSingleComplexCholesky, and IDenseComplexCholesky. They expose Size, copied L, ConditionIndicator, Solve, and SolveWithInfo. Real factors satisfy A=L*L^T; complex factors satisfy A=L*L^H.

Factorisation is O(n³) time/O(n²) factor storage; a solve with r right-hand sides is O(n²r).

QR and rank-revealing QR reference

The handles are IDenseSingleQR, IDenseDoubleQR, IDenseSingleComplexQR, and IDenseComplexQR. FactorQR and FactorPivotedQR accept m x n with m>=n.

Compact Q has shape m x n; R is n x n. For CPQR, A*P=Q*R and Permutation[j] is the source column now in factor column j. Q, R, and Permutation accessors return copies. Individual Householder column signs/phases are unspecified; reconstruction and orthogonality/unitarity are the contract.

The factor exposes Rows, Cols, NumericalRank, Tolerance, ConditionIndicator, IsColumnPivoted, SolveLeastSquares, and SolveLeastSquaresWithInfo. CPQR's rank-deficient solve is a deterministic basic solution, not a minimum-norm promise. Factorisation is O(mn²); compact factor/result storage is O(mn+n²).

The implementation is Householder QR, not classical Gram-Schmidt. CPQR chooses the largest remaining trailing column norm deterministically. See Golub and Van Loan, *Matrix Computations*, 4th ed., sections 5.2 and 5.4.

SVD and minimum-norm reference

The handles are IDenseSingleSVD, IDenseDoubleSVD, IDenseSingleComplexSVD, and IDenseComplexSVD. FactorSVD accepts tall, square, wide, singleton, and empty matrices. With p=min(m,n), U is m x p, SingularValues has length p, and V is n x p:

A=U*diag(SingularValues)*V^H.

Singular values are nonnegative and descending. Singular-vector signs/phases and bases inside repeated singular subspaces are unspecified. The factor exposes rank, tolerance, condition indicator, convergence sweep count, SolveMinimumNorm, and SolveMinimumNormWithInfo.

The deterministic one-sided Jacobi algorithm operates directly in the declared scalar type. Its leading cost is O(sweeps*m*n²), with O(mn+n²) private working storage during factorisation. See Drmač and Veselić, “New fast and accurate Jacobi SVD algorithm. I,” *SIAM J. Matrix Anal. Appl.* 29(4), 2008.

Symmetric and Hermitian eigensystem reference

FactorSymmetricEigen returns IDenseSingleSymmetricEigen or IDenseDoubleSymmetricEigen. FactorHermitianEigen returns IDenseSingleComplexHermitianEigen or IDenseComplexHermitianEigen.

Each factor exposes Size, copied ascending Eigenvalues, copied normalized column Eigenvectors, Sweeps, and Converged. Every returned factor has Converged=True; failure within the sweep limit raises instead of returning a partial result. Vector signs/phases and bases within repeated or clustered eigenspaces are unspecified. Judge results by unitary/orthogonal structure and A*v=lambda*v.

The deterministic cyclic Jacobi method costs O(sweeps*n³) and uses O(n²) factor/working storage. See Demmel, *Applied Numerical Linear Algebra*, section 5.3.

2.2 development: Hessenberg and Schur factorization

The unreleased 2.2 development unit adds real and complex overloads of ReduceHessenberg. The real factor exposes copied Q and H matrices satisfying Q^T*A*Q = H, where Q is orthogonal; the complex factor satisfies Q^H*A*Q = H, where Q is unitary. Both H factors are upper Hessenberg. The routines accept finite square matrices, including empty and singleton matrices; nil, nonsquare, non-finite, or unrepresentable reflector norms raise EDenseMatrixError. They use Householder similarity transformations, following LAPACK's documented DGEHRD and ZGEHRD contracts.

FactorRealSchur(A: IDenseDoubleMatrix) applies real shifted QR iteration and returns copied Q and T factors satisfying A = Q*T*Q^T. Q is orthogonal; T is upper quasi-triangular with isolated 1x1 real blocks and standardized 2x2 complex-conjugate blocks. The blocks are not ordered. Iterations reports the number of Francis double-shift steps. The optional MaxIterations argument sets the iteration limit; zero selects 100*max(1,n). A negative limit, invalid matrix, arithmetic overflow, or failure to converge raises EDenseMatrixError without returning a partial factor. This API provides the Schur factorization; it does not return an ordered eigenvalue list or eigenvectors.

The matrix relation, real Schur block structure, and accumulated Schur vectors follow LAPACK's documented DHSEQR contract. Single-precision Schur remains deferred. See the 2.2 design record, the real Hessenberg example, the complex Hessenberg example, and real Schur example.

FactorComplexSchur(A: IDenseComplexMatrix) returns unitary Q and upper triangular T satisfying A = Q*T*Q^H. Iterations reports the shifted-QR steps. A zero MaxIterations selects 100*max(1,n); a negative limit, invalid input, non-finite factor, or exhausted budget raises EDenseMatrixError. Because each explicit QR step costs up to O(n^3) and the iteration budget is O(n), worst-case work is O(n^4); this implementation targets moderate dense matrices. The relation and triangular form follow LAPACK's ZHSEQR contract; no LAPACK runtime is required. See the complex Schur example.

FactorRealEigen(A, Ordering, MaxIterations) builds on the real Schur factor and returns every eigenvalue as a TComplex plus its normalized right eigenvector as a column of RightEigenvectors. A real matrix can have complex conjugate pairs; those values are adjacent with the positive-imaginary member first. The default reoSchurOrder follows the Schur blocks. reoRealPart and reoMagnitude sort ascending, with ties retaining Schur order. Residuals contains the normalized backward residual for each pair, and Iterations reports the Schur steps. Eigenvectors are unit length but are not promised to be orthogonal or well-conditioned. Invalid input, arithmetic failure, or iteration exhaustion raises EDenseMatrixError without returning a partial result. The real eigenpair and conjugate-pair conventions follow LAPACK's DGEEV contract. See the runnable real nonsymmetric eigensystem example.

FactorRealGeneralizedSchur(A, B) and FactorComplexGeneralizedSchur(A, B) reduce a regular matrix pencil A - λ B. Their copied factors satisfy A = Q*S*Z^T and B = Q*T*Z^T for real input, or use conjugate transpose for complex input. The real T is upper triangular and S is upper quasi-triangular; complex S and T are upper triangular. The corresponding Alpha and Beta arrays keep each value in homogeneous form: λ = alpha/beta when beta is nonzero. A zero beta represents an eigenvalue at infinity, which can occur when B is singular. Avoid forming the quotient when its magnitude is not representable.

FactorRealGeneralizedEigen(A, B) and FactorComplexGeneralizedEigen(A, B) also return normalized right eigenvectors and normalized backward residuals for beta*A*v = alpha*B*v. Real input returns complex vectors for conjugate pairs. Accessors return defensive copies; invalid or non-finite input, indeterminate values, arithmetic failure, or iteration exhaustion raise EDenseMatrixError without a partial factor.

The implementation selects a shift gamma so A+gamma*B can be solved against B, then applies real or complex Schur reduction to the transformed matrix. It never forms B^-1*A, and keeps the projective eigenvalues as (alpha,beta) pairs. Shift selection and Schur convergence are bounded; ill-conditioned shifted solves can fail. Iterations counts the Schur stage iterations, and MaxIterations=0 selects 100*max(1,n). The factor forms and homogeneous value convention follow LAPACK's DGGES, ZGGES, and DGGEV conventions. See the generalized eigenproblem example.

Empty, non-finite, and degenerate behavior

• 0 x 0 square, QR, SVD, and eigensystem factors are valid.
• Zero right-hand-side columns produce correctly shaped empty results.
• Rank-zero matrices are valid CPQR/SVD factors; unpivoted QR solve rejects them, while SVD returns the zero minimum-norm solution.
• Repeated and clustered eigenvalues may rotate their eigenvector basis.
• NaN and infinity are rejected in every coefficient and right-hand side.
• Symmetry/Hermitian and rank decisions are relative to matrix scale and scalar precision; definiteness has its own Cholesky threshold.

The full ownership, tolerance, pivot, shape, convergence, compatibility, and algorithm-provenance decisions are fixed in the 1.6 design record.