8.3.6. Linear algebra

ulab provides a focused, microcontroller-friendly subset of numpy.linalg. Everything you need for small-matrix problems – camera calibration, IMU sensor fusion, simple kinematics – is here.

The functions are accessed as np.linalg.*:

from ulab import numpy as np

inv  = np.linalg.inv(matrix)
det  = np.linalg.det(matrix)
evs  = np.linalg.eig(matrix)

The available functions are:

  • np.linalg.cholesky(a) – Cholesky decomposition of a symmetric positive-definite matrix.

  • np.linalg.det(a) – determinant.

  • np.linalg.eig(a) – eigenvalues and eigenvectors of a real symmetric matrix.

  • np.linalg.inv(a) – matrix inverse.

  • np.linalg.norm(a) – 2-norm of a vector or matrix.

  • np.linalg.qr(a, mode='reduced') – QR decomposition.

Plus the matrix-product functions, available at the np top level:

  • np.dot(a, b) – matrix / vector product;

  • np.cross(a, b) – 3-D vector cross product;

  • np.trace(a) – sum of diagonal elements.

8.3.6.1. inv: matrix inverse

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

The inverse is computed by Gaussian elimination, so inv raises ValueError when the matrix is singular (a diagonal entry becomes zero during elimination). The cost in RAM is roughly twice the size of the input, and in time is roughly proportional to the number of entries (a 2x2 ~ 65 us, 4x4 ~ 105 us, 8x8 ~ 300 us on STM32-class hardware).

If you need to solve a linear system, do not invert and multiply – use the dedicated solvers in scipy on the OpenMV Cam (scipy.linalg.solve_triangular and scipy.linalg.cho_solve) which are both faster and numerically better behaved.

8.3.6.2. det: determinant

a = np.array([[1, 2], [3, 4]], dtype=np.uint8)
print(np.linalg.det(a))            # -2.0

The result is always a float, regardless of the input dtype. The implementation re-uses inv’s elimination, so the runtime is essentially the same.

8.3.6.3. cholesky: Cholesky decomposition

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)
# L is lower triangular; L @ L.T == a

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

8.3.6.4. eig: eigenvalues and eigenvectors

eig works only on real symmetric matrices. (Non-symmetric matrices raise ValueError.) The function returns a 2-tuple of (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)

print('eigenvalues:\n', x)
print('\neigenvectors (one per row):\n', y)

A few notes:

  • Eigenvalues are not necessarily returned in any particular order. If you need them sorted, call np.sort on the result (and apply the same permutation to the eigenvectors).

  • An eigenvector is determined only up to a non-zero scalar, so the signs of individual eigenvectors may differ from those returned by CPython numpy. This is harmless.

  • The implementation uses Givens rotations, with no closed-form estimate of run time – it iterates to convergence.

8.3.6.5. norm: vector or matrix 2-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.4162...
np.linalg.norm(m)        # 16.881...

8.3.6.6. qr: QR decomposition

qr factors a rectangular matrix A (shape (M, N)) into an orthonormal Q and an upper-triangular R such that A == Q @ R. The default mode='reduced' gives compact matrices; mode='complete' gives full (M, M) and (M, N) factors:

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

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

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

8.3.6.7. dot, cross, trace

These three live at the top of numpy, not under linalg:

from ulab import numpy as np

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

8.3.6.8. A worked example: solving a small system

To solve a 4x4 system A x = b using inversion:

from ulab import numpy as np

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

For triangular A or systems with a Cholesky factorisation already in hand, prefer scipy.linalg.solve_triangular and scipy.linalg.cho_solve – see scipy on the OpenMV Cam.

8.3.6.9. API reference