Epsilon.LinearAlgebra 1.5.0

Epsilon.LinearAlgebra

Linear algebra for Epsilon. Matrices and vectors of symbolic expressions, exact rationals or doubles: exact determinants, inverses and solutions of symbolic systems, exact elimination for rationals, LU, QR, Cholesky, eigenvalue and singular value decompositions for doubles, least squares for every entry type, and vector geometry. For games and graphics, fixed-size vectors and matrices of float, double or any number type, quaternions, camera matrices and conversions to System.Numerics. Symbolic formulas become fast code: compiled delegates or standalone C#, with Jacobians for kinematics.

var m = Matrix<Expr>.Parse("[[a, b], [c, d]]");

m.Determinant().Print();   // a * d - b * c
m.Adjugate().Print();      // [[d, -b], [-c, a]]

Installation

dotnet add package Epsilon.LinearAlgebra

Requires .NET 10. The core package Epsilon comes with it; the namespaces are Epsilon.LinearAlgebra and Epsilon.Core.

Quick start

using Epsilon.Core;
using Epsilon.LinearAlgebra;

// Symbols
var m = Matrix<Expr>.Parse("[[a, b], [c, d]]");
m.Determinant().Print();                        // a * d - b * c
m.Trace().Print();                              // a + d

// Exact numbers: the entries are rationals, never rounded
var r = Matrix<Expr>.Parse("[[1, 2], [3, 4]]");
r.Inverse().Print();                            // [[-2, 1], [3/2, -1/2]]
r.CharacteristicPolynomial("t").Print();        // t^2 - 5t - 2

// Doubles
var n = Matrix<double>.FromRows([2, 1], [4, 4]);
n.Determinant();                                // 4
n.Solve([4, 12]);                               // [1, 2]

// Vectors
Vector<double> u = [3, 4];
u.Norm();                                       // 5
Vector.Create<double>(1, 0, 0).Cross([0, 1, 0]);   // [0, 0, 1]

Features

Matrices and vectors

Matrix<T> is an immutable matrix of any entry type; Vector<T> is a column vector, a matrix with one column, and converts to one implicitly. Build them from rows, from a function of the indices, from text, or with a collection expression:

Matrix<double>.FromRows([1, 2, 3], [4, 5, 6]).Transpose();   // [[1, 4], [2, 5], [3, 6]]
Matrix<double>.Create(2, 3, (i, j) => 10 * i + j);            // [[0, 1, 2], [10, 11, 12]]
Vector<double> v = [1, 2, 3];
v.Dot([4, 5, 6]);                                            // 32
v[..2];                                                      // [1, 2]

Ranges give submatrices, and FromBlocks assembles a matrix from blocks row by row:

var g = Matrix<double>.FromRows([1, 2], [3, 4]);
g[1.., ..];                                         // [[3, 4]]
Matrix<double>.FromBlocks([g, Vector.Create(5.0, 6.0)]);   // [[1, 2, 5], [3, 4, 6]]

The arithmetic is defined for double, Rational and Expr entries: +, -, the matrix product, scalars, Pow, Hadamard, Identity, Zero, FromDiagonal, Trace. With expressions every entry of a result is simplified:

(2 * m).Print();                                 // [[2a, 2b], [2c, 2d]]
(m * Vector<Expr>.Parse("[x, y]")).Print();     // [a * x + b * y, c * x + d * y]
Matrix<Expr>.Parse("[[1, x], [0, 1]]").Pow(10).Print();   // [[1, 10x], [0, 1]]

Vector geometry

Norm, Normalize, Distance, Cross, Angle, ProjectOnto, Reflect, Lerp and Outer for vectors of doubles and of expressions. The numeric ones hold up at the extremes of floating point; the symbolic ones are exact:

Vector.Create<double>(1, -1).Reflect([0, 1]);          // [1, 1]
Vector.Create(3e300, 4e300).Norm();                    // 5E+300, not infinity
Vector.Create<double>(1, 0).Angle([1, 1e-10]);         // 1E-10, where acos gives 0
Vector<Expr>.Parse("[1, 0]").Angle(Vector<Expr>.Parse("[1, 1]")).Print();   // π / 4

Symbolic results that hold wherever they are defined

The determinant and the characteristic polynomial come from Berkowitz's algorithm, which never divides, so no entry has to be assumed nonzero. The inverse is the adjugate divided by the determinant: it is undefined exactly where the matrix is singular, and nowhere else.

var inverse = Matrix<Expr>.Parse("[[x, 1], [1, x]]").Inverse();
inverse.Print();             // [[x / (x^2 - 1), -1 / (x^2 - 1)], [-1 / (x^2 - 1), x / (x^2 - 1)]]
inverse.Substitute("x", 0).Simplify().Print();  // [[0, 1], [1, 0]]

Elimination with x as the first pivot would have divided by x and lost the point x = 0. Terms cancel exactly, functions included:

Matrix<Expr>.Parse("[[x + 1, x - 1], [x - 1, x - 3]]").Determinant().Print();    // -4
Matrix<Expr>.Parse("[[cos(x), -sin(x)], [sin(x), cos(x)]]").Inverse().Print();
                                                 // [[cos(x), sin(x)], [-sin(x), cos(x)]]

Numeric results that know when they can't be trusted

For doubles, LU decomposition with partial pivoting gives the determinant, the inverse, solutions and the rank. A pivot that is no larger than the rounding error of the terms it was computed from counts as zero, so a matrix that is singular stays singular:

var s = Matrix<double>.FromRows([1, 2, 3], [4, 5, 6], [7, 8, 9]);
s.Determinant();   // 0, not the 6.7e-16 that plain elimination gives
s.Rank();          // 2
s.Inverse();       // InvalidOperationException: The matrix is singular: its rank is 2, not 3, to working precision.

A badly scaled matrix such as [[1e-20, 0], [0, 1]] is still regular: nothing cancels there.

To compare results, IsApproximately measures the difference against the norm of the whole matrix, so the 6.1E-17 that cos(pi/2) leaves where a rotation has 0 doesn't matter:

var turn = Matrix<double>.FromRows([Math.Cos(Math.PI / 2), -1], [1, Math.Cos(Math.PI / 2)]);
turn.IsApproximately(Matrix<double>.FromRows([0, -1], [1, 0]));   // true

Least squares and decompositions

LeastSquares fits data and solves systems of any shape: by QR for doubles, exactly for rationals, and as a formula for expressions. QR, Cholesky, the eigenvalues of symmetric matrices and the singular values come with the pseudo-inverse and the condition number:

// The line y = c0 + c1 x through (0, 1), (1, 3), (2, 4) and (3, 4)
Matrix<double>.FromRows([1, 0], [1, 1], [1, 2], [1, 3]).LeastSquares([1, 3, 4, 4]);   // [1.5, 1]

var (values, vectors) = Matrix<double>.FromRows([2, 1], [1, 2]).SymmetricEigen();     // values: [1, 3]
Matrix<double>.FromRows([3, 0], [4, 5]).Svd().S;          // [6.708203932499369, 2.2360679774997894]
Matrix<double>.FromRows([1, 2], [2, 4]).PseudoInverse();  // [[0.04, 0.08], [0.08, 0.16]]
Matrix<double>.FromRows([1, 2], [3, 4]).ConditionNumber();   // 14.933034373659256

Matrix3x3<T> has the eigenvalues too, without allocating, for the principal axes of an inertia tensor or a cloud of points. The Jacobi methods behind the eigenvalues and the singular values keep small values accurate, and every test of rank or definiteness is to working precision and doesn't depend on the scale of the columns.

Exact rational matrices

Matrix<Rational> computes with exact fractions by Gaussian elimination: no rounding, no tolerance, and much faster than the same numbers as expressions. Even the ill-conditioned Hilbert matrix inverts exactly:

var hilbert = Matrix<Rational>.Create(8, 8, (i, j) => new Rational(1, i + j + 1));
hilbert.Inverse()[7, 7];                                       // 176679360
hilbert * hilbert.Inverse() == Matrix<Rational>.Identity(8);   // true
Matrix<Rational>.FromRows([1, 2], [3, 4]).Inverse();           // [[-2, 1], [3/2, -1/2]]

Fixed-size vectors and matrices for games and geometry

Vector2<T> to Vector4<T> and Matrix2x2<T> to Matrix4x4<T> are immutable value types that never allocate. Vectors are columns, so a transform reads from right to left: scale, then turn, then move:

using Vec3 = Epsilon.LinearAlgebra.Vector3<double>;
using Mat4 = Epsilon.LinearAlgebra.Matrix4x4<double>;

var model = Mat4.Translation((10, 0, 5)) * Mat4.RotationY(Math.PI / 2) * Mat4.Scale(2);
Vec3 tip = model.TransformPoint((1, 0, 0));          // [10, 0, 3]
model.Inverse().TransformPoint(tip);                 // [1, 0, 0], up to rounding
new Vec3(1, -1, 0).Reflect(Vec3.UnitY);              // [1, 1, 0]
Matrix4x4<float>.Translation((1, 2, 3)).ToSystemNumerics();   // the same map in System.Numerics

Quaternions turn and interpolate, and the camera matrices take the scene to the screen:

var turn = Quaternion<double>.RotationY(Math.PI / 2);
Quaternion<double>.Slerp(Quaternion<double>.Identity, turn, 0.5);   // an eighth of a turn
var view = Mat4.LookAt((0, 2, 5), (0, 0, 0), Vec3.UnitY);
var projection = Mat4.Perspective(Math.PI / 3, 16.0 / 9, 0.1, 100);
(projection * view).TransformPoint((0, 0, 0)).Z;    // 0.9824128790613561: the depth of the origin

The entries can be float, double, Half, int or any other number type with the arithmetic operators; norms, rotations and the inverse are there for the floating-point types. Euler angles take an explicit axis order, Compose and TryDecompose build and split translation, rotation and scale, and the projections come right- and left-handed, with the depth range of Direct3D, OpenGL or reverse-Z. The inverse knows when a matrix is singular, as the LU decomposition does: the determinant 1.7E-17 that rounding leaves for [[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]] counts as 0. Transforming many points is about as fast as with System.Numerics.

Entry-wise operations, text and LaTeX

Simplify, Expand, Substitute, Differentiate, Evaluate and EvaluateComplex work on every entry. Print() writes text that Parse reads back; ToLatex() writes a bmatrix:

var rotation = Matrix<Expr>.Parse("[[cos(x), -sin(x)], [sin(x), cos(x)]]");
rotation.Differentiate("x").Print();   // [[-sin(x), -cos(x)], [cos(x), -sin(x)]]
m.Evaluate(("a", 1), ("b", 2), ("c", 3), ("d", 4));    // [[1, 2], [3, 4]]
Vector<Expr>.Parse("[x^2, 1/2]").Simplify().ToLatex();
                                       // \begin{bmatrix} x^{2} \\ \frac{1}{2} \end{bmatrix}

From formulas to fast code

Derive a formula once with symbols, then evaluate it as fast as handwritten code. Compile turns a matrix or vector of expressions into a delegate with exactly the values of Evaluate, Jacobian gives the matrix of partial derivatives, and ToCSharp writes the formula as a C# method for code without Epsilon. Symbolic rotations and translations build chains of joints, such as an arm for inverse kinematics:

// Each link turns by the angle of its joint, then reaches along its length
Matrix<Expr> Link(string angle, int length) =>
    Matrix<Expr>.Rotation(new Variable(angle)).ToHomogeneous() * Matrix<Expr>.Translation([length, 0]);

Vector<Expr> tip = (Link("a", 3) * Link("b", 2) * Link("c", 1)).TransformPoint([0, 0]);
string[] joints = ["a", "b", "c"];
var position = tip.Compile(joints);
var jacobian = tip.Jacobian(joints).Compile(joints);
position([0, 0, 0]);           // [6, 0]: stretched out along x
jacobian([0, 0, 0]).Row(1);    // [6, 3, 1]: how far up each joint moves the tip per radian
File.WriteAllText("ArmTip.cs", tip.ToCSharp("ArmTip", joints));   // public static void ArmTip(double a, ...

A subexpression shared by several entries is computed once: a compiled 3 x 3 rotation by three angles takes about 60 nanoseconds instead of 2 microseconds. The guide solves the inverse kinematics of this arm in a few more lines.

Good to know

• System.Numerics has a Vector<T> of its own, so Vector<double> is ambiguous in a file that uses both namespaces. Name the type you need with an alias, using RealVector = Epsilon.LinearAlgebra.Vector<double>;, or write System.Numerics.BigInteger in full instead of importing the namespace.
• Vector3<T> and Matrix4x4<T> don't clash with Vector3 and Matrix4x4 of System.Numerics, which have no type parameter. System.Numerics multiplies row vectors, so its matrices are the transposes of these; ToSystemNumerics() and FromSystemNumerics() transpose for you.
• Like ExprParser.Parse, Parse doesn't simplify: [[1/2, x + x]] keeps 1 / 2 and x + x until Simplify().
• The operators on Expr matrices simplify in the default Generic mode, so an entry x/x becomes 1.

Documentation

The guides explain every feature in detail, the API reference lists every public type and member, and the changelog lists the changes. The expressions themselves are documented with the core.

License

MIT, copyright (c) 2026 Max Zakharov - see LICENSE.