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.

// Delphi — same function names, free functions (not class methods)
y := Mul(A, x);           // Matrix × vector → Vector
y := Mul(x, A);           // Vector × matrix → Vector
C := Mul(A, B);           // Matrix × matrix → Matrix
C := Mul(x, y);           // Outer product → Matrix
R := Mul(A, B, C, D);     // Chain product

With TVec/TMtx (method calls, zero allocation)

// Matrix × matrix: C.Mul(A, B)
C.Mul(A, B);                                     // C := A * B
C.Mul(A, B, opTranspose);                        // C := A^T * B
C.Mul(A, B, mtSymmetric, mtGeneral);             // exploit symmetry

// Matrix × vector: result.TensorProd(Mtx, Vec) or result.TensorProd(Vec, Mtx)
y.TensorProd(A, x);                              // y := A * x
y.TensorProd(x, A);                              // y := x * A

// Outer product: C.TensorProd(colVec, rowVec)
C.TensorProd(x, y);                              // C := x * y^T

// Element-wise matrix multiply (Hadamard product):
C.MulElem(A, B);                                 // C := 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 := A * B * C
R.Mul(A, B, C, D);                           // R := A * B * C * D

// Arbitrary-length chain — pass the factors as an open array literal.
// No temporary variable, no SetLength.
R.MulArray([M0, M1, M2, M3, M4]);            // R := M0 * M1 * M2 * M3 * M4

// Optional transpArray pre-transposes individual factors without
// allocating a transposed temporary. `1` means transpose, `0` means
// use as-is; the literal position matches the factor position.
R.MulArray([M0, M1, M2, M3, M4], [0, 1, 0, 0, 1]);
// equivalent to 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:

// Two-arg: operate on self
self.MulDiagLeft (v);         // self := diag(v) * self   (row-scale)
self.MulDiagRight(v);         // self := self * diag(v)   (column-scale)

// Three-arg: explicit source matrix, store into self
self.MulDiagLeft (v, M);      // self := diag(v) * M
self.MulDiagRight(M, v);      // self := M * diag(v)

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:

// Compute R := A * diag(w) * B  without allocating a diagonal matrix.
tmp.MulDiagRight(A, w);       // tmp := A * diag(w)  (column-scale A)
R.Mul(tmp, B);                // R := tmp * B

// Or fold into a longer chain:
R.MulArray([A, B, C, D]);     // R := A*B*C*D
R.MulDiagRight(w);             // R := R * diag(w)   (post-scale columns)

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):

// Matrix times matrix
S.MulLeft (X, Y);             // Y := S * X       (X, Y: TMtx)
S.MulRight(X, Y);             // Y := X * S       (X, Y: TMtx)

// Matrix times vector
S.MulLeft (x, y);             // y := S * x       (x, y: TVec)
S.MulRight(x, y);             // y := x * S       (x, y: TVec)

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)

x := LDivide(A, b);      // Solve A * x = b
x := Divide(b, A);       // Solve x * A = b

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, mtSymmetric);             // exploit symmetry

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

// SVD-based (minimum-norm solution)
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 then ...                     // 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, mtSymmetric);                    // symmetric → real eigenvalues
A.Eig(VL, D, VR);                         // with eigenvectors

MtxExpr convenience wrappers

eigenvalues := Eig(A);                    // → Vector of eigenvalues
singularvals := SVD(A);                   // → Vector of singular values
ok := Cholesky(A);                        // in-place, returns boolean

Matrix Inverse

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

// TMtx method (in-place)
A.Inv;                                    // A is replaced by A^(-1)
A.Inv(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

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.

In Delphi, Vector/Matrix DO benefit from the object cache due to compiler-managed lifetime. Operator expressions are efficient in Delphi.

// Operators are efficient in Delphi (object cache handles temporaries)
c := a + b;           // temporary from cache, returned automatically
d := c * e;           // same — no GC pressure

// TVec methods are still faster for hot loops (zero allocation)
c.Add(a, b);
c.Mul(e);

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)`