LSQSolver

日本語版

A lightweight least-squares solver for dense matrices in .NET. It can be used for linear and polynomial regression, curve fitting with user-defined basis functions, parameter estimation from redundant observations, and inverse problems with fewer observations than unknowns. If your model is linear in its unknown parameters, it can usually be written as a least-squares problem and solved with LSQSolver.

Supports:

• Overdetermined systems
• Underdetermined systems
• Rank-deficient systems
• Minimum-norm solutions
• Multiple right-hand sides

without SVD and without external dependencies.

Related projects: LSQSolver · LSQSolver.Complex · LSQSolver.MathNet

---

Features

• Column-Pivoted QR Decomposition (CPQR)
• Automatic numerical rank detection
• Minimum-norm solutions
• Rank-deficient problem support
• Parallelized factorization using Parallel.For
• Optional low-allocation overwrite mode
• Optional storage of QR intermediates (R, Qᵀb, and pivot information)
• Column-major array API for direct integration with numerical code
• Input validation and non-throwing status reporting for numerical failures
• No external dependencies

---

Why LSQSolver?

Many numerical libraries solve least-squares problems through SVD-based pseudo-inverses.

While SVD is extremely robust, it can be more expensive than necessary when the goal is simply to obtain a least-squares solution.

LSQSolver is based on:

1. Column-Pivoted QR Decomposition (CPQR)
2. Numerical rank detection
3. Cholesky-based minimum-norm reconstruction

This approach supports rank-deficient and underdetermined problems without requiring a full singular value decomposition.

---

Installation

dotnet add package LSQSolver

Usage

Import module

using LSQSolver;
using static LSQSolver.LSQSolver;

Construct a matrix

double[][] data =
{
    new double[] { 1, 2 },
    new double[] { 3, 4 },
    new double[] { 5, 6 }
};

var A = new MatrixObject(data);

Solve a least-squares problem

double[] b = { 7, 8, 9 };
var result = Solve(A, b); // or LSQSolver.LSQSolver.Solve(A, b);

For example, rows of A may represent observations and columns may represent regression terms, basis functions, or unknown model parameters.

Solve Method

var result = Solve(
    A,
    b,
    overwrite: true,
    store_intermediates: false,
    rank_tolerance: 2.22044604925032e-16,
    check_finite: true);

The same options are available for all overloads. The main overloads are:

Solve(double[] columnMajorMatrix, int rows, int cols, double[] b, ...);
Solve(double[] columnMajorMatrix, int rows, int cols, double[] B, int rhs_count, ...);
Solve(MatrixObject A, double[] b, ...);
Solve(MatrixObject A, double[] B, int rhs_count, ...);
Solve(MatrixObject A, MatrixObject B, ...);

Parameters

• A / columnMajorMatrix: Coefficient matrix. The array overload expects a rows-by-cols matrix in column-major order, where (i, j) is stored at j * rows + i.
• b / B: One right-hand side, or a rows-by-rhs_count column-major matrix of right-hand sides.
• rhs_count: Number of right-hand sides stored in B.
• overwrite: If true, the supplied matrix array and right-hand side storage are used as work arrays and overwritten.
• store_intermediates: If true, stores R, Qᵀb, and pivot information.
• rank_tolerance: Relative tolerance used for numerical rank detection.
• check_finite: If true, checks whether A and b contain NaN or Infinity.

Return value: LSQSolverResult

Solve() returns an LSQSolverResult object.

The result object contains the computed solution, solver status, basic diagnostics, and optional QR intermediate data.

Before using Solution, check Status.

Property	Type	Description
Status	LSQSolverStatus	Status code of the solve operation. Check this value before using Solution.
Solution	double[]	Computed solution in column-major order (Cols-by-RHSCount). For rank-deficient or underdetermined problems, each column is the minimum-norm solution when the solve succeeds.
SolutionMatrix	MatrixObject?	Solution represented as a Cols-by-RHSCount matrix. This is especially convenient for multiple right-hand sides.
Rows	int	Number of rows of the input matrix A.
Cols	int	Number of columns of the input matrix A.
RHSCount	int	Number of right-hand sides solved simultaneously.
Rank	int	Estimated numerical rank detected during the column-pivoted QR factorization.
ResidualNorm	double	Euclidean residual norm for one right-hand side, or Frobenius norm ||AX - B||_F for multiple right-hand sides.
RArray	double[]?	Column-major storage containing the CPQR result. Only its upper-triangular/trapezoidal part represents R; entries below it retain internal Householder data. Stored only when store_intermediates is true.
R	MatrixObject?	Obsolete compatibility view of the same storage as RArray. New code should use RArray. The same warning about entries below the upper-triangular part applies.
Qtb	double[]?	Transformed right-hand side Qᵀb or matrix QᵀB, stored column-major. This is stored only when store_intermediates is true; otherwise it is null.
Pivot	int[]?	Column pivot information. Pivot[j] gives the original column index of the j-th pivoted column. This is stored only when store_intermediates is true; otherwise it is null.
Tag	object?	Optional tag field reserved for additional metadata.

Public methods

Method	Description
ToString()	Returns the complete result in a readable form.
ToString(bool omit, int display_row_count = 10, int display_col_count = 10)	Returns a readable result preview. If omit is true, matrix output is limited to the requested row and column counts.

Example

var result = Solve(A, b, store_intermediates: true);

if (result.Status != LSQSolverStatus.Success)
{
    Console.WriteLine(result.Status);
    return;
}

double[] x = result.Solution;

Console.WriteLine($"Rank: {result.Rank}");
Console.WriteLine($"Residual norm: {result.ResidualNorm}");
Console.WriteLine(result.ToString(omit: true));

Multiple right-hand sides

When the same coefficient matrix is used with several data vectors, place them in the columns of B and solve them together:

var B = new MatrixObject(new double[][]
{
    new double[] { 7, 1 },
    new double[] { 8, 0 },
    new double[] { 9, 1 }
});

var result = Solve(A, B);
if (result.Status == LSQSolverStatus.Success)
{
    MatrixObject X = result.SolutionMatrix!; // A.Cols x B.Cols
}

The factorization of A is shared across all right-hand sides. The array overload can be used to avoid conversion when A and B are already stored in column-major arrays.

Column-major array API

double[] a =
{
    1, 3, 5, // first column
    2, 4, 6  // second column
};
double[] b = { 7, 8, 9 };

var result = Solve(a, rows: 3, cols: 2, b: b);

Documentation

Theory

theory.md

Topics:

• Least-squares formulation
• Rank-revealing QR decomposition
• Numerical rank detection
• Minimum-norm solutions
• Cholesky-based reconstruction

Examples

polynomial-fit.md

• Polynomial curve fitting
• Practical fitting examples

gravity-inversion.md

• Gravity anomaly inversion
• Underdetermined least-squares problems
• Minimum-norm reconstruction

Performance

Preliminary benchmark results are available in:

performance.md

The benchmark compares LSQSolver with GNU Octave for dense square least-squares problems under both full-rank and rank-deficient conditions.

The reported timings are medians of 10 runs and were measured on the following machine:

Item	Value
Model	MacBook Pro
Chip	Apple M1 Pro
CPU Cores	8 total: 6 performance cores and 2 efficiency cores
Memory	16 GB

In this benchmark, LSQSolver showed the following results:

• For full-rank dense matrices, LSQSolver was faster than Octave QR factorization for n >= 50 in this benchmark.
• At n = 2000, LSQSolver took 564.4 ms for the full-rank case, while Octave QR factorization took 1483.0 ms.
• For rank-deficient matrices, LSQSolver became faster than Octave QR factorization for larger sizes. At n = 2000, LSQSolver took 764.7 ms, while Octave QR factorization took 2102.2 ms.
• Compared with Octave pinv, LSQSolver was significantly faster for large matrices. At n = 2000, LSQSolver was about 12.8x faster for the rank-deficient case and about 27.6x faster for the full-rank case.

These results are preliminary and depend on hardware, runtime, compiler settings, BLAS/LAPACK configuration, and benchmark implementation details.

Development Direction

LSQSolver is intended to remain a focused, lightweight dense least-squares solver rather than grow into a general-purpose numerical framework. Future extensions will be evaluated according to whether they:

• cover a recurring least-squares use case that is difficult or error-prone to assemble correctly on the caller side;
• avoid unnecessary copies or repeated factorizations by being implemented inside the solver; or
• provide a clear numerical or algorithmic advantage over a straightforward transformation of the input problem.

Weighted and regularized least-squares support remain possible directions, but will not be added solely to expand the API. Backward compatibility, failure safety, predictable storage, and maintainable numerical kernels take priority over feature count.

Related Projects

Project	Description
LSQSolver	The core rank-aware least-squares solver for real-valued dense problems.
LSQSolver.Complex	Complex-valued least-squares support built on LSQSolver.
LSQSolver.MathNet	MathNet.Numerics integration for real and complex least-squares problems.

License

MIT License