Linear Algebra vs Element-wise Operations

The * and / operators on Vector/Matrix are always element-wise — they are not matrix multiplication or linear system solving. This is the #1 source of confusion.

Several method names are overridden on TMtx/Matrix so that on a two-dimensional container they perform the real linear-algebra operation instead of the element-wise one their base-class signature implies. The key overrides are:

• TMtx.Mul(A, B) with two TMtx arguments → matrix multiplication (self := A * B), not element-wise.
• TMtx.Inv → matrix inverse (self := self^-1), not per-element 1/x. On TVec/Vector the inherited Inv is per-element.

So whether .Mul/.Inv is element-wise or linear-algebra depends on the receiver type: element-wise on TVec/Vector, linear-algebra on TMtx/Matrix — the polymorphism is deliberate.

The Critical Distinction

Operation	Element-wise form	Linear-algebra form
Multiply — * operator	c := a * b on Vector/Matrix (always per-element)	— (operators never do linear algebra)
Multiply — method	c.Mul(a, b) on TVec/Vector (per-element)	C.Mul(A, B) on TMtx/Matrix (matrix product) — overridden, see above
Divide	c := a / b on Vector/Matrix (per-element)	MtxExpr.Divide(b, A) — solves A^T · x = b
Back-divide	—	MtxExpr.LDivide(A, b) — solves A · x = b
Inverse	v.Inv on TVec/Vector (per-element 1/x)	M.Inv on TMtx/Matrix (matrix inverse) — overridden

Two API Layers

Vector/Matrix — value types with operator overloading, return new objects:

c := a + b;     // allocates result from object cache
c := a * b;     // ELEMENT-WISE multiply (not matrix product!)
c := a * 2.0;   // scalar multiply

TVec/TMtx — pool-managed types, write result into calling object (self-as-destination):

// On TVec (lowercase v, a, b are TVec — vectors):
v.Add(a, b);       // v := a + b       (no allocation)
v.Mul(a, b);       // v := a .* b      (element-wise — the only meaning on 1D)
v.Mul(2.0);        // v := v * 2.0     (scalar multiply)
v.Inv;             // v := 1 ./ v      (element-wise 1/x)

// On TMtx (uppercase M, A, B are TMtx — matrices) the same method
// names are OVERRIDDEN to perform real linear algebra:
M.Add(A, B);       // M := A + B       (element-wise; addition is the same)
M.Mul(A, B);       // M := A * B       (matrix product — not element-wise)
M.MulElem(A, B);   // M := A .* B      (explicit element-wise product on TMtx)
M.Inv;             // M := M^-1        (matrix inverse — not per-element)

The receiver type decides the meaning: on a 1D container (TVec/Vector) every .Mul/.Inv is element-wise because there is nothing else it can be. On a 2D container (TMtx/Matrix) .Mul(A, B) and .Inv mean real linear algebra. Use .MulElem on TMtx when you actually want the Hadamard (per-element) product.

Matrix Multiplication

With Vector/Matrix (via MtxExpr static functions)

MtxExpr functions accept TVec/TMtx inputs and return Vector/Matrix value types. These are the primary way to do linear algebra with expression-friendly syntax.

// C#
using Dew.Math.Units;

Vector y = MtxExpr.Mul(A, x);           // A * x  (matrix × vector)
Vector y = MtxExpr.Mul(x, A);           // x * A  (vector × matrix)
Matrix C = MtxExpr.Mul(A, B);           // A * B  (matrix × matrix)
Matrix C = MtxExpr.Mul(x, y);           // x * y^T (outer product)

// With transpose/conjugate
Matrix C = MtxExpr.Mul(A, B, TMtxOperation.opTranspose);  // A^T * B

// Chain: A * B * C * D
Matrix R = MtxExpr.Mul(A, B, C, D);
// or: MtxExpr.Mul(new TMtx[] { A, B, C, D });

With TVec/TMtx (method calls, zero allocation)

// C# — identical syntax (TMtx/TVec are the same API)
C.Mul(A, B);
y.TensorProd(A, x);
y.TensorProd(x, A);
C.TensorProd(x, y);
C.MulElem(A, B);

Chain multiplication — one call, any number of matrices

TMtx.Mul has overloads for three- and four-factor chains, and TMtx.MulArray accepts an arbitrary TMtxArray so the whole chain is computed in a single call. These are preferred over manual chaining (tmp.Mul(A,B); C.Mul(tmp,D);) because the implementation picks the most efficient evaluation order and avoids redundant temporaries.

// Three and four matrices in one call
R.Mul(A, B, C);
R.Mul(A, B, C, D);

// Arbitrary chain via TMtx[]
TMtx[] chain = new TMtx[] { M0, M1, M2, M3, M4 };
R.MulArray(chain);                           // R := M0*M1*M2*M3*M4

int[] transp = { 0, 1, 0, 0, 1 };            // transpose M1 and M4
R.MulArray(chain, transp);                   // R := M0 * M1^T * M2 * M3 * M4^T

Chain factors must all be matrices. The element type of TMtxArray is TMtx, and the fixed-arity Mul(A, B, C[, D]) overloads take TMtx arguments. There is no heterogeneous [TMtx, TVec, TMtx, ...] overload. For vector weights and vectors-as-operands see the next two subsections.

Folding a vector weight into a chain — `MulDiagLeft` / `MulDiagRight`

When one of the factors in a chain is a diagonal matrix whose diagonal is known as a TVec, you don't need to materialise the full diagonal matrix. TMtx exposes four overloads that multiply in a vector treated as the main diagonal — equivalent to row or column scaling:

From the summaries: "DiagMtx contains the values of the main diagonal of the diagonal matrix. … This operation is the same as calling TMtx.ScaleRows" (left form) or "TMtx.ScaleCols" (right form). The vector length must match self.Cols for MulDiagLeft and self.Rows for MulDiagRight.

Example — A · diag(w) · B in two calls, no diagonal matrix built:

This is strictly faster than materialising diag(w) as a square TMtx and including it in a MulArray — the scaling step is O(n²) instead of the O(n³) a full matrix multiply would cost.

Sparse matrix multiplication — `TSparseMtx.MulLeft` / `MulRight`

TSparseMtx (sparse container) does not share TMtx's Mul/MulArray API. Instead it exposes MulLeft / MulRight, each with a matrix-matrix and a matrix-vector overload. These are the only dense-API methods that take a TVec operand directly as a vector (not as a diagonal):

Summaries from the DB: "Multiply the sparse matrix from left" / "Multiply sparse matrix from left with vector X and place the result in vector Y" — and the symmetrical right-hand forms. Use these when your matrix is stored as a TSparseMtx; for dense matrix-vector products use TVec.TensorProd(Mtx, Vec) as shown above.

Note: TSparseMtx has no chain overload — MulLeft/MulRight are always binary. For a sparse chain you apply them one factor at a time.

Solving Linear Systems (A·x = b)

With MtxExpr (returns Vector/Matrix)

// Solve A * x = b  (left-divide, like MATLAB's A\b)
Vector x = MtxExpr.LDivide(A, b);         // uses LQR internally
Matrix X = MtxExpr.LDivide(A, B);         // multi-column RHS

// Solve x * A = b  (right-divide, like MATLAB's b/A)
Vector x = MtxExpr.Divide(b, A);

With TMtx (method calls, more control)

// LU-based (fastest for square, well-conditioned)
A.LUSolve(b, x);                          // A * x = b, general
A.LUSolve(b, x, TMtxType.mtSymmetric);    // exploit symmetry

// QR-based (robust, handles rank deficiency)
A.LQRSolve(b, x);                         // full rank
int rank = A.LQRSolve(b, x, 1e-9);        // rank-deficient, returns effective rank

// SVD-based (minimum-norm solution)
int rank = A.SVDSolve(b, x, s, 1e-9);     // s receives singular values

Decompositions

With TMtx (in-place, reuse objects)

// Cholesky (symmetric positive definite)
if (A.Cholesky()) { ... }                 // A is replaced by factor

// LU decomposition
A.LU(LU_result, pivots);

// QR / LQ decomposition
A.LQR(L, Q, R);                           // economy-size by default
A.LQR(L, Q, R, P);                        // with column pivoting

// SVD
A.SVD(U, S, V);                           // full SVD
A.SVD(S);                                 // singular values only

// Eigenvalues
A.Eig(D);                                 // eigenvalues in D
A.Eig(D, TMtxType.mtSymmetric);           // symmetric → real eigenvalues
A.Eig(VL, D, VR);                         // with eigenvectors

MtxExpr convenience wrappers

Vector eigenvalues = MtxExpr.Eig(A);              // eigenvalues
Vector singularvals = MtxExpr.SVD(A);              // singular values
bool ok = MtxExpr.Cholesky(A);                     // in-place, returns boolean

Matrix Inverse

// MtxExpr (returns Matrix value type)
Matrix B = MtxExpr.Inverse(A);                           // general inverse
Matrix B = MtxExpr.Inverse(A, TMtxType.mtSymmetric);     // exploit symmetry

// TMtx method (in-place)
A.Inv();                                   // A is replaced by A^(-1)
A.Inv(TMtxType.mtSymmetric);              // symmetric inverse

Important: Prefer LDivide(A, b) over Inverse(A) * b. Solving is faster and more numerically stable than computing the explicit inverse.

Performance: Vector/Matrix vs TVec/TMtx

In C#, Vector, Matrix, VectorInt, and MatrixInt do not benefit from the object cache. C# language design prevents user-defined memory allocation control for value types. Every operator expression allocates new managed objects. Use TVec/TMtx with CreateIt/FreeIt for best performance.

lifetime. Operator expressions are efficient in Delphi.

// SLOW — Vector allocates on every operation
Vector c = a + b;           // allocates new Vector
Vector d = c * e;           // allocates another

// FAST — TVec reuses pooled memory
TVec c = TVec.CreateIt(a.Length);
try {
    c.Add(a, b);            // no allocation
    c.Mul(e);               // in-place, no allocation
} finally {
    TVec.FreeIt(ref c);
}

// Linear algebra — same trade-off:
Vector y = MtxExpr.Mul(A, x);   // convenient but allocates
y.TensorProd(A, x);              // fast, writes into pre-allocated y

Quick Reference

The * and / operators on Vector/Matrix are always element-wise — never matrix multiplication. Use the method calls or MtxExpr free functions for linear algebra.

Task	Vector/Matrix operator	MtxExpr / free function	TMtx/TVec method
Element-wise add / sub	c := a + b, a - b	—	c.Add(a, b), c.Sub(a, b)
Element-wise multiply	c := a * b	—	c.Mul(a, b) on TVec; C.MulElem(A, B) on TMtx
Element-wise divide	c := a / b	—	c.DivideVec(a, b) (per-element); see also DivideBy(N) for N / self
Element-wise 1/x	—	—	v.Inv on TVec/Vector
Matrix × matrix	— (operators are element-wise!)	Mul(A, B) → Matrix	C.Mul(A, B) on TMtx
Matrix × vector	—	Mul(A, x) → Vector	y.TensorProd(A, x)
Vector × matrix	—	Mul(x, A) → Vector	y.TensorProd(x, A)
Outer product x · yᵀ	—	Mul(x, y) → Matrix	C.TensorProd(x, y)
Row-scale (diag(w) · A)	—	—	A.MulDiagLeft(w) ≡ ScaleRows
Column-scale (A · diag(w))	—	—	A.MulDiagRight(w) ≡ ScaleCols
Sparse mat × mat / mat × vec	—	—	S.MulLeft(X, Y) / MulRight(X, Y) on TSparseMtx
Chain matrix product	—	Mul(A, B, C[, D]) → Matrix	R.Mul(A, B, C[, D]), R.MulArray([...])
Solve A · x = b	—	LDivide(A, b) → Vector	A.LUSolve(b, x) or A.LQRSolve(b, x)
Solve x · A = b	—	Divide(b, A) → Vector	—
Matrix inverse	—	Inverse(A) → Matrix	A.Inv (overridden on TMtx)
Eigenvalues	—	Eig(A) → Vector	A.Eig(D)
SVD	—	SVD(A) → Vector	A.SVD(U, S, V)