8.13. Linear algebra¶

Linear algebra on a camera is small-matrix work: a 3x3 rotation that fuses an IMU sample into the world frame, a calibration matrix that rectifies a lens, a Kalman filter’s state-covariance update, a polynomial fit whose normal equations come out as a tiny linear solve. The numpy.linalg and scipy.linalg submodules cover exactly that scale.

The functions live in two modules. Matrix-product operations are at the numpy top level; decompositions and the matrix inverse are under numpy.linalg; the dedicated linear-system solvers are under scipy.linalg. A 2-by-2 matrix multiply, for example:

from ulab import numpy as np

A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
#        [43.0, 50.0]])

8.13.1. What is available¶

dot() – matrix or vector product.
cross() – 3-D vector cross product.
trace() – sum of the diagonal.
inv() – matrix inverse.
det() – determinant.
cholesky() – Cholesky decomposition (symmetric positive-definite input).
eig() – eigenvalues and eigenvectors of a real symmetric matrix.
norm() – 2-norm of a vector or matrix.
qr() – QR decomposition with mode='reduced' (default) or mode='complete'.
solve_triangular() – solve A @ x = b when A is triangular.
cho_solve() – solve A @ x = b given a Cholesky factor of A.

8.13.2. dot, cross, trace¶

dot() is how matrix multiplication is written on the camera. The @ operator that desktop numpy provides for the same job is not implemented on ndarray, so every matrix-vector, matrix-matrix, and vector dot product goes through the dot() call:

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.dot(a, b)             # 32.0 (scalar product)
np.cross(a, b)           # array([-3.0, 6.0, -3.0])

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a)             # matrix-vector product
np.dot(m, m)             # matrix-matrix product
np.trace(m)              # 1 + 5 + 9 = 15.0

The result of dot() is always of dtype float. cross() is the cross product of two 3-vectors, trace() the sum of the main diagonal of a square matrix.

8.13.3. inv and det¶

m = np.array([[1, 2, 3, 4],
              [4, 5, 6, 4],
              [7, 9, 9, 4],
              [3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))

The inverse is computed by Gauss-Jordan elimination, so inv() raises ValueError when the matrix is singular (a diagonal entry becomes zero during elimination). The RAM cost is roughly twice the size of the input.

The determinant re-uses the same elimination – runtime is essentially the same as the inverse.

When the application’s goal is to solve a linear system, do not invert and multiply – prefer the dedicated solvers below. Both are faster and better-behaved numerically.

8.13.4. cholesky¶

For a symmetric positive-definite matrix A, cholesky() returns a lower-triangular L such that A = L @ L.T:

a = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]])
L = np.linalg.cholesky(a)

If the input is not positive definite or not symmetric, ValueError is raised.

The Cholesky factor is half the work of an LU factorisation and is the right starting point for any problem whose matrix is known to be symmetric positive definite (covariance updates, normal equations from a least-squares fit).

8.13.5. eig¶

eig() works only on real symmetric matrices. Non-symmetric matrices raise ValueError. It returns a 2-tuple (eigenvalues, eigenvectors):

a = np.array([[1, 2, 1, 4],
              [2, 5, 3, 5],
              [1, 3, 6, 1],
              [4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)

Notes:

The eigenvalues come back in no particular order. Apply sort() (and the same permutation to the eigenvectors via argsort()) when a sorted order matters.
An eigenvector is unique only up to a non-zero scalar, so the sign of individual eigenvectors is not uniquely defined. Two correct runs may produce vectors of opposite sign; this is harmless.

8.13.6. norm¶

The Euclidean (Frobenius) norm of a vector or matrix:

v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

np.linalg.norm(v)        # 7.416...
np.linalg.norm(m)        # 16.881...

The optional axis= keyword takes the norm along a single axis instead of over the whole array.

8.13.7. qr¶

qr() factors a rectangular matrix A (shape (M, N)) into an orthonormal Q and an upper-triangular R such that A == Q @ R:

A = np.arange(6).reshape((3, 2))

q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)

q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)

The decomposition is implemented via successive Givens rotations. The right choice for least-squares problems where the matrix is not symmetric.

8.13.8. Solving systems¶

The two dedicated solvers under ulab.scipy.linalg are both faster and more accurate than np.dot(np.linalg.inv(A), b):

solve_triangular(a, b, lower=False)() – solve a @ x = b assuming a is triangular:

A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

cho_solve(L, b)() – given a Cholesky factor L, solve A @ x = b where A = L @ L.T:
```
L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)
```

Reach for these instead of inverting whenever the structure of A allows – they save elimination work and the explicit inverse.

8.13.9. Solving a small linear system¶

A = np.array([[3, 0, 1, 1],
              [0, 1, 0, 2],
              [1, 0, 1, 1],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])

x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x))         # should equal b

The same problem can be expressed by factoring A and calling the appropriate solver – faster and more accurate when A has the right structure.

For the complete argument-level reference, see numpy.linalg — Linear algebra routines and scipy.linalg — Linear algebra routines.