# How to solve simultaneous equations faster with Compile?

I have large 6x6 matrix Uwhich is a multiplication of 15 rotational matrix. All of the elements are Sin\[theta] and Cos\[Theta]. I have to solve simultaneous equations using the matrix elements. Yesterday i ran the code for 14 hours with core i3 machine at windows 7. Still no output. I know that Complie can run code quicker but i am not an expert on this. So i need help. here is the code.I tried but failed to use complie here. how to use complie or how can i run the code faster?

Solve[U[[2,4]]==0.12 && U[[1,4]]==0.11 && U[[2,5]]==0.17 && U[[1,5]]==0.11 &&
U[[2,6]]==0.14 && U[[1,6]]==0.11 &&
U[[3,1]]==0    && U[[3,2]]==0    && U[[3,3]]==0    && U[[3,4]]==0    &&
U[[3,5]]==0    && U[[3,6]]==0,
{\[Theta]12, \[Theta]13, \[Theta]14, \[Theta]15, \[Theta]16, \[Theta]23, \[Theta]24,
\[Theta]25, \[Theta]26, \[Theta]34, \[Theta]35, \[Theta]36, \[Theta]45, \[Theta]46,
\[Theta]56}
]

-
Generally,you should try to provide self-contained examples so that the code can be run, to maximize our chances to provide reasonable answers. Here it means providing your matrix U. –  Leonid Shifrin Sep 28 '12 at 10:38
The short answer is that you can't expect Compile to help you here. Solve is not compilable and it very unlikely that you could write your own linear solver that does a better job than Solve. If you have a 16x16 linear system and decimal approximations are sufficient, then Solve should be able to solve it nearly instantaneously. The information you have as written is really not quite sufficient to advise further, though. –  Mark McClure Sep 28 '12 at 12:02

There exists no solution. A product of rotation matrices is invertible, but the last 6 constraints introduce an entire row of zeros in U, making its rank of 5 or less, whence it must not be invertible.

Here, then, is optimal code: :-)

Solve[False, {}]


(0 seconds).

Well, maybe that wasn't so constructive. Let's see how we can find a product of rotation matrices in $n$ dimensions with specified components. This method minimizes the $L^2$ norm of the differences at the desired locations.

Begin by specifying the places and the values to attain at those places:

places = {{2, 4}, {1, 4}, {2, 5}, {1, 5}, {2, 6}, {1, 6},
{3, 1}, {3, 2}, {3, 3}, {3, 4}, {3, 5}, {3, 6}};
targets = {.12, .11, .17, .11, .14, .11, 0, 0, 0, 0, 0, 0};


The objective function is the sum of squares of differences at those places:

objective[u_] := With[{x = u~Extract~places - targets}, x.x];


(This clear formulation is thanks to Mr.Wizard.)

Rotation matrices in higher dimensions, for rotations from one basis vector to another, look like the usual rotation matrices in two dimensions, with ones elsewhere on the diagonal:

rotationMatrix[t_, n_Integer, i_Integer, j_Integer] /; 1 <= i <= n && 1 <= j <= n && i != j :=
Module[{u = IdentityMatrix[n], s = Sin[t], c = Cos[t]},
u[[{i, j}, {i, j}]] = {{c, -s}, {s, c}}; u
]


The variables to parameterize a generic rotation correspond to these "basis vector" rotations, which when multiplied give an arbitrary rotation. The variables, vars, can be created at the same time the matrix product u is generated, cached via Sow, and collected via Reap:

With[{n = 6},
vars = Flatten[Last[Reap[
u = Dot @@ (Flatten[
Table[rotationMatrix[Sow[Unique["\[Theta]"]], n, i, j],
{i, 1, n - 1}, {j, i + 1, n}], 1])
]]]
];


Let's give it a whirl:

solution = NMinimize[objective[u], vars]


(0.69 seconds).

In this case, the best we can do is reduce the objective function to $1.$, not $0.$ (which would indicate a good solution). We might characterize what we have found as a least-squares fit to the targets. For the record, here it is:

Chop[u /. Last[solution]] // MatrixForm


-
I'd replace Outer[0 &, Range[n], Range[n]] with ConstantArray[0, {n, n}] or even Table[0, {n}, {n}]. –  Guess who it is. Sep 28 '12 at 16:18
@J.M. Thanks--I knew there was a simple way, but I just wanted to get the job done and this method worked fine. I'm sure others (and maybe you again) will have additional improvements to suggest. (E.g., you might find a nicer way to express objective.) If you like, just edit the post to make them; if not, I'll wait a bit and then edit it all at once. –  whuber Sep 28 '12 at 16:20
It seems as though objective can be written: objective[u_] := With[{x = u ~Extract~ places - targets}, x.x] or objective = #.# &[u ~Extract~ places - targets] & –  Mr.Wizard Sep 28 '12 at 23:07