Solving linear systems

Sooner or later every robotics problem turns into the same question: find the unknown vector x such that A * x = b. Node voltages in a circuit, forces in a frame, joint angles that reach a point, coefficients that calibrate a sensor: all of them are linear systems. Octave has a single operator for it, the backslash, and learning when to trust its answer is worth more than learning the operator itself.

From a circuit to Ax = b

Take a ladder of three nodes. Node 1 has a 2 mS path to ground, node 2 has 1 mS, node 3 has 1 mS. Nodes 1 and 2 are joined by 1 mS (a 1 kilo-ohm resistor), and so are nodes 2 and 3. A source pushes 5 mA into node 1 and 1 mA into node 3. Kirchhoff's current law says the currents leaving a node add up to the current injected there, so for each node we write one equation, using conductance times voltage difference:

node 1: 2*v1 + 1*(v1 - v2) = 5 which is 3*v1 - v2 = 5

node 2: 1*v2 + 1*(v2 - v1) + 1*(v2 - v3) = 0 which is -v1 + 3*v2 - v3 = 0

node 3: 1*v3 + 1*(v3 - v2) = 1 which is -v2 + 2*v3 = 1

Collect the coefficients and you get the conductance matrix G. Its diagonal holds the total conductance touching a node, and every off-diagonal entry is minus the conductance of the link between two nodes. Using millisiemens and milliamps keeps the numbers small and returns volts directly. The same pattern describes a chain of springs in a structure: a stiffness matrix times displacements equals forces.

The backslash operator

G = [3 -1 0; -1 3 -1; 0 -1 2];   % conductance matrix, mS
i_src = [5; 0; 1];               % injected currents, mA
v = G \ i_src;
printf("node voltages (V): %.3f %.3f %.3f\n", v);
printf("det(G) = %.3f\n", det(G));
printf("rank(G) = %d\n", rank(G));
disp(norm(G*v - i_src) < 1e-12)

node voltages (V): 2.000 1.000 1.000
det(G) = 13.000
rank(G) = 3
1

Check the answer by substitution, as you always should: node 1 gives 3*2 - 1 = 5, node 2 gives -2 + 3 - 1 = 0, node 3 gives -1 + 2 = 1. The last line does the same check numerically: the residual G*v - i_src is essentially zero, so disp prints 1 for true.

You read G \ i_src as "G divided into i_src from the left". Octave factors G (an LU factorisation with row pivoting) and solves by substitution. It never builds the inverse. People new to the language write inv(G) * i_src, which does more work and is less accurate, because forming the inverse amplifies rounding error that the factorisation avoids. Remember the rule: if you see inv multiplied by a vector, you almost certainly want a backslash.

Determinant and rank

The determinant of G is 13. It is not zero, so the system has exactly one solution. The rank is the number of independent equations, and rank 3 for three unknowns confirms it.

When rank falls below the number of unknowns, the matrix is singular and the answer is not unique. A classic physical cause is a floating node: a part of the circuit with no path to ground. Its voltage can be shifted up or down without changing any current, so no single answer exists. Octave tells you with warning: matrix singular to machine precision. A simpler example is S = [1 2; 2 4], where row 2 is twice row 1 so rank(S) returns 1: the second equation carries no new information.

Do not use det as a singularity test on real data. The determinant depends on units and size: ten equations with a scale of 0.1 have determinant 1e-10 and are perfectly healthy. Use rank or, better, the condition number below.

Conditioning

A system can be non-singular and still be treacherous. Take two nearly parallel lines: x1 + x2 = 2 and x1 + 1.0001*x2 = b2. The condition number cond(A) measures how much a relative error in b can be magnified in x. Roughly, you lose log10(cond(A)) of the 16 or so decimal digits a double carries.

A = [1 1; 1 1.0001];
printf("cond(A) = %.0f\n", cond(A));
for b2 = [2.0001 2.00015]
  x = A \ [2; b2];
  printf("b2 = %.5f -> x = %.4f %.4f\n", b2, x);
end

cond(A) = 40002
b2 = 2.00010 -> x = 1.0000 1.0000
b2 = 2.00015 -> x = 0.5000 1.5000

Changing one measurement by 0.00005, a relative change of about 0.0025 percent, moves the solution from (1, 1) to (0.5, 1.5). That is a 50 percent change in the answer. A sensor with noise at that level would make the result meaningless, and nothing is wrong with the algorithm: the geometry of the problem is to blame. Robots meet this constantly. Two range beacons seen from almost the same direction, calibration points that lie almost on a line, and an arm stretched out near a singular pose all give nearly parallel equations.

The same maths in NumPy

numpy_version.py

ResetRun

import numpy as np # Node conductances in mS and injected currents in mA, so volts come out directly G = np.array([[ 3.0, -1.0, 0.0], [-1.0, 3.0, -1.0], [ 0.0, -1.0, 2.0]]) i_src = np.array([5.0, 0.0, 1.0]) v = np.linalg.solve(G, i_src) print("node voltages (V):", " ".join("%.3f" % x for x in v)) print("det(G) = %.3f" % np.linalg.det(G)) print("rank(G) =", np.linalg.matrix_rank(G)) print("residual small:", np.linalg.norm(G @ v - i_src) < 1e-12) S = np.array([[1.0, 2.0], [2.0, 4.0]]) print("rank of singular S =", np.linalg.matrix_rank(S)) A = np.array([[1.0, 1.0], [1.0, 1.0001]]) print("cond(A) = %.0f" % np.linalg.cond(A)) for b2 in (2.0001, 2.00015): x = np.linalg.solve(A, np.array([2.0, b2])) print("b2 = %.5f -> x = %.4f %.4f" % (b2, x[0], x[1]))

np.linalg.solve(G, i_src) is the NumPy counterpart of the backslash: same factorisation, same result. The runner prints the same node voltages, determinant and condition number as Octave. Try changing 1.0001 to 1.001 in the matrix: cond(A) drops by a factor of ten, to about 4000.