Regularization of singular linear system

Question

I have

$$\begin{pmatrix} 0&B_3&-B_2 \\ -B_3&0&B_1 \\ B_2&-B_1&0 \end{pmatrix}\begin{pmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{pmatrix}=\begin{pmatrix} \Delta_1 \\ \Delta_2 \\ \Delta_3 \end{pmatrix}$$

With $\Delta_1 B_1 + \Delta_2 B_2 + \Delta_3 B_3 = 0$. Because the constraint, the system has clearly solutions in the null space of the $B$ matrix for every possible set of values $\Delta$ that satisfy the constraint.

However, I'm not sure what is the best way to incorporate the constraint into LinearSolve so that it will be able to find the whole set of solutions.

I tried replacing one of the $\Delta$ in terms of the others, but no matter which $\Delta$ I choose to replace in terms of the others, I keep getting solutions where $\omega_3$ will be zero. This seems weird since nothing in the problem makes a special distinction between the last coordinate. I've noticed that RowReduce would turn the matrix into

$$\begin{pmatrix} 1&0&-B_1/B_3 \\ 0&1&-B_2/B_3 \\ 0&0&0 \end{pmatrix}$$

In any case, if I replace $\Delta_1$ with the constraint equation, I get a solution

$$\begin{pmatrix} -\Delta_2 /B_3 \\ -B_2 \Delta_2/B_1 B_3 - \Delta_3/B_1 \\ 0 \end{pmatrix} $$

Intuitively, the plane orthogonal to $B$ as a vector contains the elements of $\Delta$ and also a set of solutions to the inhomogeneous system of equations, which can be extended with NullSpace in order to construct the whole set of solutions. But I wonder why these solutions are constrained to the $\omega_3 = 0$ subspace in the first place

Code that manifests the issue:

LinearSolve[{{0,B3, -B2},{-B3,0,B1},{B2, -B1,0}},{-d2*(B2/B1)-d3*(B3/B1),d2,d3}]

Also

LinearSolve[{{0,B3, -B2},{-B3,0,B1},{B2, -B1,0}},{d1,-d1*(B1/B2)-d3*(B3/B2),d3}]

Answer 1

Let's create some example data:

d = {d1, d2, d3};
A = {{B1, B2, B3}};
d = d - PseudoInverse[A].A.d // ComplexExpand // Simplify;
B = {{0, B3, -B2}, {-B3, 0, B1}, {B2, -B1, 0}};
A.d // Simplify

(* {0} *)

Now, we augment the system matrix B and the right hand side d:

BB = ArrayFlatten[{{B, - Transpose[A]}, {A, 0}}];
dd = Join[d, {0}];
x = LinearSolve[BB, dd][[1 ;; 3]]
(* {(-B3 d2 + B2 d3)/(B1^2 + B2^2 + B3^2), (B3 d1 - B1 d3)/(B1^2 + B2^2 + B3^2), (-B2 d1 + B1 d2)/(B1^2 + B2^2 + B3^2)} *)

Let's test the result:

B.x - d // Simplify
A.x // Simplify

(* {0,0,0} *)
(* {0} *)

With the augmentation, we enforced that $B_1 \, \omega_1 + B_2 \, \omega_2 + B_3 \, \omega_3 = 0$ and we added the orthogonal complement of the image of B to the image of BB.

Answer 2

Your problem appears ill posed. Call the matrix mat, and define

mat = {{0, b3, -b2}, {-b3, 0, b1}, {b2, -b1, 0}};
b = {b1, b2, b3};
d = {d1, d2, d3};

You have two conditions: one is that b.d==0. The other is that d==mat.w, where w is the unknown. Substituting d into the first equation shows that b.mat.w==0. But this is true for all w, since b.mat is identically zero!

b.mat

{0, 0, 0}