Skip to content

deepbox/linalg

Solvers

Dense, triangular, banded, least-squares and sparse solvers, and the Sylvester and Lyapunov equations.

  • Linear systems

Overview

  • solve(a, b) uses LU with partial pivoting and throws DataValidationError when a is singular. lstsq(a, b) returns { x, residuals, rank, s } and gives the minimum-norm solution for rank-deficient or underdetermined systems. Unlike NumPy, residuals is always filled in.
  • solveTriangular(a, b, lower = true) treats a as lower triangular by default, where SciPy defaults to upper. Pass false for an upper matrix. The fourth argument takes { trans, unitDiagonal }.
  • sparseSolve and sparseCholeskySolve accept a CSRMatrix from deepbox/ndarray or the result of denseToCSR. sparseSolve expands the matrix into a dense workspace, so memory is O(N^2) and time is O(N^3) in the worst case. It is not a sparse direct solver. solveBanded solves banded systems directly and is the better choice for tridiagonal matrices.
  • sylvester(a, b, c) solves AX + XB = C and throws when A and -B share an eigenvalue. lyapunov(a, q) solves AX + XAᵀ = Q.

Types

3

CSRMatrix

type

Compressed Sparse Row (CSR) matrix representation.

Stores a sparse matrix using three arrays: - values: Non-zero values (NNZ elements) - colIndices: Column index for each non-zero value (NNZ elements) - rowPointers: Start index in values/colIndices for each row (N+1 elements)

For row i, the non-zero entries are at indices rowPointers[i]..rowPointers[i+1]-1. Repeated (row, column) entries are summed, as in scipy.sparse.

This record is not the CSRMatrix class of deepbox/ndarray, but the solvers accept both (see SparseMatrixInput).

  • readonly n: number
  • readonly values: Float64Array
  • readonly colIndices: Int32Array
  • readonly rowPointers: Int32Array

SolveTriangularOptions

type

Options for solveTriangular.

  • readonly trans?: boolean

    Solve A^T x = b instead of A x = b (default: false).

  • readonly unitDiagonal?: boolean

    Treat the diagonal of A as all ones and do not read it (default: false).

SparseMatrixInput

type

A sparse matrix accepted by sparseSolve and sparseCholeskySolve: either the plain CSRMatrix record made by denseToCSR, or the CSRMatrix class of deepbox/ndarray (which must be square).

type SparseMatrixInput = CSRMatrix | NdarrayCSRMatrix;

Functions

9

denseToCSR

function denseToCSR(a: Tensor): CSRMatrix;

Convert a dense tensor to CSR sparse format.

Parameters

a
Tensor
Input 2D tensor (must be square)

Returns

CSRMatrixCSR representation

Throws

  • ShapeErrorIf input is not a square 2-D matrix
  • DTypeErrorIf input has string or complex dtype
  • DataValidationErrorIf input contains non-finite values (NaN, Infinity)
Example
import { tensor } from 'deepbox/ndarray';import { denseToCSR } from 'deepbox/linalg'; const A = tensor([[4, 1, 0], [1, 3, 1], [0, 1, 2]]);const csr = denseToCSR(A);

lstsq

function lstsq(a: Tensor, b: Tensor, rcond?: number): { readonly x: Tensor; readonly residuals: Tensor; readonly rank: number; readonly s: Tensor; };

Least squares solution to A * x = b.

Finds x that minimizes ||A*x - b||^2 (Euclidean norm). Works for overdetermined (M > N), underdetermined (M < N), and square systems. For rank-deficient A the minimum-norm solution is returned.

Algorithm: SVD-based least squares

Parameters:

Parameters

a
Tensor
Coefficient matrix of shape (M, N)
b
Tensor
Target values of shape (M,) or (M, K)
rcondoptional
number
Cutoff for small singular values: singular values <= rcond * s_max are treated as zero (default: float64 machine epsilon * max(M,N), as in NumPy, which always solves in double precision) Returns: Object with: - x: Least squares solution of shape (N,) or (N, K) - residuals: Sum of squared residuals ||b - A*x||^2, shape (1,) for 1-D b or (K,) for 2-D b. Unlike NumPy, it is always filled in, also for rank-deficient or underdetermined systems. - rank: Effective rank of A - s: Singular values of A Requirements: - A can be any M x N matrix - First dimension of A must match first dimension of b

Returns

{ readonly x: Tensor; readonly residuals: Tensor; readonly rank: number; readonly s: Tensor; }

Throws

  • ShapeErrorIf A is not 2D matrix
  • ShapeErrorIf b is not 1D or 2D tensor
  • ShapeErrorIf dimensions don't match
  • DTypeErrorIf input has string or complex dtype
  • InvalidParameterErrorIf rcond is negative or non-finite
  • DataValidationErrorIf input contains non-finite values (NaN, Infinity)
Example
import { lstsq } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; // Overdetermined system (more equations than unknowns)const A = tensor([[1, 1], [1, 2], [1, 3]]);const b = tensor([2, 3, 5]);const result = lstsq(A, b); console.log(result.x);  // Best fit solutionconsole.log(result.residuals);  // Residual error

lyapunov

function lyapunov(a: Tensor, q: Tensor): Tensor;

Solve the continuous Lyapunov equation AX + XAᵀ = Q.

This is a special case of the Sylvester equation with B = Aᵀ.

Parameters

a
Tensor
Square matrix A of shape (N, N)
q
Tensor
Symmetric matrix Q of shape (N, N)

Returns

TensorSolution matrix X of shape (N, N)

Example
import { lyapunov } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[-1, 0], [0, -2]]);const Q = tensor([[2, 0], [0, 8]]);const X = lyapunov(A, Q);// A*X + X*A^T ≈ Q

solve

function solve(a: Tensor, b: Tensor): Tensor;

Solve linear system A * x = b.

Finds vector x that satisfies the equation A * x = b.

Algorithm: LU decomposition with partial pivoting

Parameters:

Parameters

a
Tensor
Coefficient matrix of shape (N, N)
b
Tensor
Right-hand side of shape (N,) or (N, K) for multiple RHS Returns: Solution x of same shape as b Requirements: - A must be square matrix - A must be non-singular (invertible) - Number of rows in A must equal length of b Numerical singularity: as in NumPy, only an exactly zero pivot is reported as singular. A matrix that is singular up to rounding error returns a solution with very large entries instead of throwing. Check cond(a) or matrixRank(a) first when the input may be rank deficient, or use lstsq().

Returns

Tensor

Throws

  • ShapeErrorIf A is not square
  • ShapeErrorIf dimensions don't match
  • DTypeErrorIf input has string or complex dtype
  • DataValidationErrorIf A is singular
  • DataValidationErrorIf input contains non-finite values (NaN, Infinity)
Example
import { solve } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[3, 1], [1, 2]]);const b = tensor([9, 8]);const x = solve(A, b); // Verify: A @ x ≈ bconsole.log(x);  // [2, 3]

solveBanded

function solveBanded(luBands: readonly [ number, number ], ab: Tensor, b: Tensor): Tensor;

Solve a banded linear system A x = b.

Details

The banded matrix A is specified by its lower and upper bandwidths and the band data in compact banded storage (same convention as scipy.linalg.solve_banded).

Band storage format: ab has shape (l + u + 1, N) where ab[u + i - j, j] = A[i, j] for max(0, j-u) <= i <= min(N-1, j+l). Entries of ab that do not correspond to a matrix element (the corners) are ignored, but must still be finite.

Tridiagonal systems that are diagonally dominant use the Thomas algorithm (no pivoting is needed for those); every other system uses banded LU with partial pivoting. The matrix is factored once, however many right-hand sides there are.

Time Complexity: O(N * (l + u)²) for the factorization plus O(N * (l + u)) per right-hand side, much faster than dense O(N³)

Parameters

luBands
readonly [ number, number ]
Tuple [l, u] where l = number of lower diagonals, u = number of upper diagonals
ab
Tensor
Band matrix in compact form, shape (l + u + 1, N)
b
Tensor
Right-hand side vector of shape (N,) or matrix of shape (N, nrhs)

Returns

TensorSolution x with same shape as b (float64)

Throws

  • InvalidParameterErrorIf l or u is not a non-negative integer
  • ShapeErrorIf ab is not 2-D with l + u + 1 rows, or b does not match N
  • DTypeErrorIf ab or b has string or complex dtype
  • DataValidationErrorIf the matrix is singular, or ab or b contain NaN or Infinity
Example
import { solveBanded } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; // Tridiagonal system: [2 -1 0; -1 2 -1; 0 -1 2] x = [1; 0; 1]// l=1, u=1, band storage://   row 0 (upper): [0, -1, -1]//   row 1 (diag):  [2,  2,  2]//   row 2 (lower): [-1, -1, 0]const ab = tensor([[0, -1, -1], [2, 2, 2], [-1, -1, 0]]);const b = tensor([1, 0, 1]);const x = solveBanded([1, 1], ab, b); // [1, 1, 1]

solveTriangular

function solveTriangular(a: Tensor, b: Tensor, lower?: boolean, options?: SolveTriangularOptions): Tensor;

Solve triangular system.

More efficient than solve() when A is already triangular. Only the selected triangle of A is used; the entries of the other triangle are ignored (they must still be finite).

Parameters:

Parameters

a
Tensor
Triangular matrix of shape (N, N)
b
Tensor
Right-hand side of shape (N,) or (N, K)
loweroptional
boolean
If true (default), A is lower triangular; if false, upper triangular. Note that scipy.linalg.solve_triangular defaults to upper.
optionsoptional
SolveTriangularOptions
trans solves with the transpose of A (which swaps the triangle that is read), unitDiagonal assumes a unit diagonal Returns: Solution x of the same shape as b (float64)

Returns

Tensor

Throws

  • ShapeErrorIf A is not square
  • ShapeErrorIf dimensions don't match
  • DTypeErrorIf input has string or complex dtype
  • DataValidationErrorIf A is singular (zero diagonal; ignored with unitDiagonal)
  • DataValidationErrorIf input contains non-finite values (NaN, Infinity)
Example
import { solveTriangular } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const L = tensor([[2, 0], [3, 4]]);  // Lower triangularconst b = tensor([6, 18]);const x = solveTriangular(L, b, true);               // [3, 2.25]const y = solveTriangular(L, b, true, { trans: true }); // solves L^T y = b

sparseCholeskySolve

function sparseCholeskySolve(input: SparseMatrixInput, b: Tensor): Tensor;

Solve a sparse symmetric positive-definite system using sparse Cholesky.

Details

Performs Cholesky decomposition A = L Lᵀ where L is lower triangular, then solves via forward/backward substitution.

Only the lower triangle of the sparse matrix is read (the upper triangle may be absent or hold anything); symmetry is not checked. L is stored in envelope (skyline) form: row i keeps the columns from its first stored entry to the diagonal, which is where Cholesky fill-in is confined. Work and memory therefore scale with the envelope, not with N².

Time Complexity: O(N * bandwidth²) for banded SPD matrices Space Complexity: O(envelope size), at most O(N²) for dense rows

Parameters

input
SparseMatrixInput
b
Tensor
Right-hand side vector of shape (N,)

Returns

TensorSolution vector x of shape (N,)

Throws

  • ShapeErrorIf b is not a vector of length N, or rowPointers does not have N + 1 entries
  • DataValidationErrorIf the CSR structure is invalid, or the matrix is not positive definite
Example
import { tensor } from 'deepbox/ndarray';import { denseToCSR, sparseCholeskySolve } from 'deepbox/linalg'; const A = tensor([[4, 2, 0], [2, 5, 1], [0, 1, 3]]);const b = tensor([1, 2, 3]);const x = sparseCholeskySolve(denseToCSR(A), b);

sparseSolve

function sparseSolve(input: SparseMatrixInput, b: Tensor): Tensor;

Solve a sparse linear system A x = b by LU factorization with partial pivoting.

Details

The CSR matrix is expanded into a dense n x n workspace and factored with the same LU routine as solve, skipping the elimination of zero entries. This is fast for moderately sized sparse systems, but memory use is O(N²) and time is O(N³) in the worst case: it is not a sparse direct solver with fill-reducing ordering. For symmetric positive-definite systems with a small bandwidth, sparseCholeskySolve stores only the envelope.

Time Complexity: O(N²) to scan for pivots plus the cost of the non-zero updates (O(N * bandwidth²) for banded matrices), O(N³) for dense ones Space Complexity: O(N²)

Parameters

input
SparseMatrixInput
b
Tensor
Right-hand side vector of shape (N,)

Returns

TensorSolution vector x of shape (N,)

Throws

  • ShapeErrorIf b is not a vector of length N, or rowPointers does not have N + 1 entries
  • DataValidationErrorIf the CSR structure is invalid (index out of range, decreasing rowPointers, non-finite values), or the matrix is singular
Example
import { tensor } from 'deepbox/ndarray';import { denseToCSR, sparseSolve } from 'deepbox/linalg'; const A = tensor([[4, 1, 0], [1, 3, 1], [0, 1, 2]]);const b = tensor([1, 2, 3]);const x = sparseSolve(denseToCSR(A), b);

sylvester

function sylvester(a: Tensor, b: Tensor, c: Tensor): Tensor;

Solve the continuous Sylvester equation AX + XB = C.

Details

Uses the Bartels–Stewart algorithm: 1. Compute real Schur decompositions A = U*Ta*U^T and B = V*Tb*V^T 2. Transform to U^T*C*V = Ta*Y + Y*Tb where Y = U^T*X*V 3. Solve the quasi-triangular system column by column 4. Recover X = U*Y*V^T

The equation has a unique solution only when A and -B share no eigenvalue. When they do, or come within rounding error of doing so (relative to the largest entry of the Schur forms), a DataValidationError is thrown instead of returning a meaningless huge solution.

Time Complexity: O(M³ + N³ + M²N) where A is M×M and B is N×N

Parameters

a
Tensor
Square matrix A of shape (M, M)
b
Tensor
Square matrix B of shape (N, N)
c
Tensor
Matrix C of shape (M, N)

Returns

TensorSolution matrix X of shape (M, N) (float64)

Throws

  • ShapeErrorIf a, b or c is not 2-D, a or b is not square, or c is not (M, N)
  • DTypeErrorIf an input has string or complex dtype
  • DataValidationErrorIf an input contains NaN or Infinity, or if A and -B share an eigenvalue (the solution is not unique)
Example
import { sylvester } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 0], [0, 2]]);const B = tensor([[3, 0], [0, 4]]);const C = tensor([[4, 5], [6, 12]]);const X = sylvester(A, B, C);// A*X + X*B ≈ C

Deprecated names

1

These older spellings still work in 1.x. New code should use the names above.

  • solve_bandedPrefer solveBanded.
linalg-solvers.ts
import {  denseToCSR,  lstsq,  lyapunov,  solve,  solveBanded,  solveTriangular,  sparseSolve,  sylvester,} from "deepbox/linalg";import { tensor } from "deepbox/ndarray"; const a = tensor([[3, 1], [1, 2]]);const b = tensor([9, 8]); console.log(solve(a, b).toString());console.log(solveTriangular(a, b, false).toString()); // upper triangularconst fit = lstsq(a, b);console.log(fit.x.toString(), fit.rank);console.log(sparseSolve(denseToCSR(a), b).toString());console.log(sylvester(a, a, tensor([[1, 0], [0, 1]])).toString());console.log(lyapunov(tensor([[-1, 0], [0, -2]]), tensor([[2, 0], [0, 8]])).toString()); // Tridiagonal system with one lower and one upper diagonal, in band storage.const bands = tensor([[0, -1, -1], [2, 2, 2], [-1, -1, 0]]);console.log(solveBanded([1, 1], bands, tensor([1, 0, 1])).toString());

Found a mistake? Open an issue.