dsh-linalg

Linear algebra toolbox for DeepSeek Harness (dsh): matrix multiply, determinant, inverse, trace, transpose, reduced row echelon form, rank, integer matrix powers, a linear system solver that classifies unique / infinitely many / no solutions, a symmetric eigen-decomposition with measured residuals, matrix factorizations (LU with partial pivoting, thin QR, Cholesky) with measured reconstruction residuals, a singular value decomposition with rank and condition number, a least-squares solver for overdetermined systems, and vector operations (dot / cross / norm / projection / angle).

Zero runtime dependencies, pure local arithmetic — no network, no processes, no eval.

Why

Language models frequently make arithmetic errors on matrix multiplication, determinants, inverses and linear systems. This plugin hands those computations to exact local code with Gaussian elimination and partial pivoting, so results are deterministic and correct to 10 decimal places.

Install

dsh plugin --profile web add github:TYEclipse/dsh-linalg

Replace web with your profile name. (Requires pnpm on your PATH.)

Tools

Tool	What it does
matrix_multiply	Multiply two matrices A·B with dimension checking
matrix_compute	One operation per call: transpose, determinant, inverse, trace, rref, rank, power (integer exponent)
matrix_decompose	One factorization per call: lu (P·A = L·U, singular flagged), qr (thin A = Q·R, R diagonal positive) or cholesky (A = L·Lᵀ, symmetric positive definite only) — each with a measured reconstruction residual
matrix_svd	Singular value decomposition A = U·Σ·Vᵀ: all singular values (descending), numerical rank, condition number σ₁/σ_rank, thin U/V, measured residuals, and an optional Moore–Penrose pseudo-inverse
matrix_eigen	Eigen-decomposition of a real symmetric matrix: all eigenvalues (descending) + unit eigenvectors + measured residual
least_squares	Minimizer of ‖A·x − b‖ for overdetermined systems, plus the minimum residual norm and max ‖Aᵀ(A·x − b)‖ (rank-deficient / underdetermined input is an error)
solve_linear	Solve Ax = b; classifies unique / infinite (particular + nullspace basis + free-variable count) / no solution
vector_ops	dot, cross (3D), norm, projection, angle (degrees)

Input matrices are plain JSON arrays of arrays, e.g. [[1, 2], [3, 4]]. Maximum dimension is 20×20 by default (configurable).

Examples

Determinant:

matrix_compute(matrix=[[2,1,1],[1,3,2],[1,0,0]], op="determinant")  →  -1

Inverse:

matrix_compute(matrix=[[4,7],[2,6]], op="inverse")
→ [[0.6,-0.7],[-0.2,0.4]]

Solve 2x+3y=8, x−y=1:

solve_linear(a=[[2,3],[1,-1]], b=[8,1])
→ unique solution: [2.2, 1.2]

Under-determined system (1 equation, 3 variables) reports the full parametric form:

solve_linear(a=[[1,1,1],[1,2,3]], b=[6,14])
→ kind=infinite, particular=[-2,8,0], nullspaceBasis=[[1,-2,1]], freeVariableCount=1

Cross product:

vector_ops(op="cross", a=[1,2,3], b=[4,5,6])  →  [-3,6,-3]

Angle between [3,4] and the x-axis:

vector_ops(op="angle", a=[3,4], b=[4,0])  →  53.1301023542  (degrees)

Rank and integer powers:

matrix_compute(matrix=[[1,2,3],[4,5,6],[7,8,9]], op="rank")  →  2
matrix_compute(matrix=[[1,1],[1,0]], op="power", exponent=20)
→ [[10946, 6765], [6765, 4181]]

Eigenvalues of a symmetric matrix (2x2 by hand would be a quadratic; here it is exact):

matrix_eigen(matrix=[[2,1],[1,2]])
→ eigenvalues (descending): 3, 1
  sum = 4, max residual = 0
  eigenvectors (unit, one per row):
    lambda1 = 3: [0.7071067812, 0.7071067812]
    lambda2 = 1: [0.7071067812, -0.7071067812]

The 4x4 tridiagonal [[2,-1,0,0],[-1,2,-1,0],[0,-1,2,-1],[0,0,-1,2]] gives 3.6180339887, 2.6180339887, 1.3819660113, 0.3819660113 — the closed form 2 − 2·cos(kπ/5).

LU factorization (row permuted, one swap):

matrix_decompose(matrix=[[2,1,1],[4,-6,0],[-2,7,2]], kind="lu")
→ p = [1, 0, 2]  (row i of P·A is row p[i] of A)
  swaps = 1, sign = -1, singular = false, det from U = -16
  L = [[1,0,0],[0.5,1,0],[-0.5,1,1]]
  U = [[4,-6,0],[0,4,1],[0,0,1]]

Thin QR (R diagonal positive, so the factorization is unique):

matrix_decompose(matrix=[[3,0],[4,5]], kind="qr")
→ Q = [[0.6,-0.8],[0.8,0.6]]   R = [[5,4],[0,3]]

Cholesky of a symmetric positive definite matrix:

matrix_decompose(matrix=[[4,12,-16],[12,37,-43],[-16,-43,98]], kind="cholesky")
→ L = [[2,0,0],[6,1,0],[-8,5,3]], det from L = 36

Least squares — fit y = a + b·x through (1,1), (2,2), (3,2):

least_squares(a=[[1,1],[1,2],[1,3]], b=[1,2,2])
→ x1 = 0.6666666667, x2 = 0.5
  min ‖A·x − b‖ = 0.4082482905, max ‖Aᵀ(A·x − b)‖ = 0

Singular values of a 3×2 matrix, with the pseudo-inverse:

matrix_svd(matrix=[[1,2],[3,4],[5,6]], pseudoInverse=true)
→ singular values (descending): 9.5255180916, 0.5143005807
  rank 2 of 2, condition number 18.5213053413
  pseudo-inverse A⁺ (2x3): [[-1.3333333333, -0.3333333333, 0.6666666667],
                            [1.0833333333, 0.3333333333, -0.4166666667]]

Numerical behaviour

• Gaussian elimination with partial pivoting (largest-magnitude pivot per column).
• All output numbers are rounded to 10 decimal places (configurable roundPlaces); values below 5e-12 are snapped to exactly 0, so -0 never appears.
• Pivot tolerance is 1e-12: entries below that are treated as zero, so numerically singular matrices are reported as singular instead of producing garbage.
• Singular inverse requests fail cleanly with an error message; inconsistent systems return kind: "none" with the reducing row equation.
• matrix_eigen runs the cyclic Jacobi method in full double precision: every eigenvector is unit length, eigenvectors are signed so the leading non-zero component is positive, and the reply carries maxResidual = the largest |A·v − λ·v| entry actually measured. Non-square or non-symmetric input, and non-convergence, are reported as errors — never a guessed spectrum.
• power uses binary exponentiation; negative exponents invert first (and fail cleanly on singular input), exponents are limited to the integer range −64…64.
• matrix_decompose reports the residual it actually measured (max |P·A − L·U|, max |A − Q·R|, max |A − L·Lᵀ|). A singular matrix still gets an LU pair, but singular: true is set and the determinant from U is 0; linearly dependent columns make QR fail (no orthonormal completion is invented); an asymmetric or indefinite matrix is rejected by Cholesky with the offending entry / pivot spelled out.
• least_squares solves through the thin QR factorization (not the normal equations, which square the condition number) and returns both the minimum residual norm and max |Aᵀ(A·x − b)| — zero proves the answer really is the minimizer. Rank-deficient or underdetermined input is an error rather than one arbitrary minimizer.
• matrix_svd runs one-sided Jacobi: pairs of columns are rotated until they are mutually orthogonal, so the column norms are the singular values and the accumulated rotations are the right singular vectors. Rank uses the standard σ > max(m,n)·ε·σ₁ tolerance, and it is a true zero rather than a small number: [[1,2,3],[4,5,6],[7,8,9],[10,11,12]] reports rank 2 with σ₃ = 0 (forming AᵀA would square the condition number and lose that zero). A rank-deficient matrix gets no condition number — the honest answer is "infinite", reported as an absent field plus a note, not as a huge number. Factors are signed so the leading non-zero component of each right singular vector is positive, and the reply carries max |A − U·Σ·Vᵀ| plus the largest |UᵀU − I| / |VᵀV − I| entry.

Configuration

- name: 'github:TYEclipse/dsh-linalg'
  config:
    maxDimension: 20   # max rows/columns (1–50)
    roundPlaces: 10    # decimals in every output number (1–15)

Development

pnpm install
pnpm build
pnpm test
pnpm lint

License

MIT