Safely produce c = a / b
for matrices a
and b
of appropriate
dimensions. This is fit for code generation.
c = smrdivide(a, b);
We can do this by taking a singular value decomposition of b
, inverting
its singular values, transposing it, and multiplying it on the left by
a
.
[u, s, v] = svd(b);
inv(b) == v * diag(diag(s).^-1) * u.'
a * inv(b) == a * v * diag(diag(s).^-1) * u.'
Of course, to be safe, we don't want to invert the singular values when they're zero, so we check for this.
Two optional tolerances can be used to determine if the singular values are big enough for inversion: an absolute tolerance and a relative tolerance. For an absolute tolerance of, e.g., 1e-15, no singular value less than 1e-15 will be inverted. For a relative tolerance of 1e-10, no singular value less than the largest singular value times this relative tolerance will be inverted. If both are provided, the maximum tolerance of the two will be used.
Clearly, this matrix doesn't really produce a / b
-- that's undefined
where b
has null space. However, the intended use of this function is
for safety in critical applications, where b
should be full rank but
for whatever reason isn't. It keeps the program from crashing with an
error about division by 0.
a | Any matrix |
---|---|
b | The matrix by which to right-divide |
abs_tol | Absolute tolerance; singular values smaller than this number won't be inverted. Can be left empty. |
rel_tol | Relative tolerance; singular values smaller than this number times the largest singular value won't be inverted. By default, 1e-15. |
c | The result of |
---|
Make one normal matrix, a
, and a matrix with null space, b
.
a = randn(3, 5);
a(:, 2) = 0;
b = randn(5);
b(:, 2) = 0;
b(2, :) = 0;
Show that the rank of b is 4, not 5 we'd need for full inversion.
rank(b)
ans =
4
Try a normal divide -- doesn't work.
c1 = a / b
[Warning: Matrix is singular to working precision.]
[> In ml2jade (line 254)
In enjaden (line 360)
In deploy_web (line 62)
]
c1 =
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
Divide manually, ignoring null space.
c2(:, [1 3 4 5]) = a(:, [1 3 4 5]) / b([1 3 4 5], [1 3 4 5])
c2 =
0.7546 0 -0.9916 -1.1112 -0.3590
-0.8776 0 0.7010 1.1744 0.4129
0.7195 0 -0.3537 -0.7209 -0.2368
Show that the safe divide gives the same answer as ignoring null space.
c3 = smrdivide(a, b)
c3 =
0.7546 -0.0000 -0.9916 -1.1112 -0.3590
-0.8776 0.0000 0.7010 1.1744 0.4129
0.7195 -0.0000 -0.3537 -0.7209 -0.2368
*kf
v1.0.3