Calculating a matrix with a specific form

125 views (last 30 days)
Jan Buchali on 7 Apr 2021
Commented: Steven Lord on 8 Apr 2021
Hello,
I am having troubles with calculating a Matrix of a specific form out of the Equation Ax = B, where A is the Matrix im looking for, B is a 3x4728 Matrix and x is also an 3x4728 Matrix. B and x are measurments.
Thats why Im using A = B*X'*inv(X*X') for the calculation. I know that A has to be in the form of
[ 0, -m(3), m(2);
m(3), 0, -m(1);
-m(2), m(1), 0].
My Problem is now that im getting an A Matrix but not in the right form.
Does anybody has an idea how to get an Matrix A in the right form?
Steven Lord on 8 Apr 2021
Let's take a step back. You've identified an approach to solving your problem that you've asked the group to help troubleshooting. But it's possible there's a more robust and/or easier approach to solve your underlying problem than the one you've identified.
Can you describe the problem you're trying to solve that led you to this particular formulation of A*x = B? Is James Tursa correct in guessing "Is this some form of small angle approximation from measurements? Are you trying to find a small angle rotation matrix?"

James Tursa on 7 Apr 2021
Edited: James Tursa on 7 Apr 2021
You could rearrange the equations, isolate m(1), m(2), and m(3) and solve for them directly. E.g., rearrange the equations to form
Xbig * m = Bbig
then solve for m elements. E.g.,
z = zeros(size(x,2),1);
Xbig = [ z x(3,:)' -x(2,:)';
-x(3,:)' z x(1,:)';
x(2,:)' -x(1,:)' z ];
Bbig = reshape(B',[],1);
m = Xbig\Bbig;
A = [ 0, -m(3), m(2);
m(3), 0, -m(1);
-m(2), m(1), 0 ];
Caveat: This was quickly done on paper ... the code above is untested.
Side question: Is this some form of small angle approximation from measurements? Are you trying to find a small angle rotation matrix? Maybe there is a better approach you can use for finding rotation matrices from measurement data (Matt J has an FEX contribution for this).
Matt J on 7 Apr 2021
James' code looks good to me. A quick test on synthetic data,
mtrue=rand(3,1);
m=mtrue;
A = [ 0, -m(3), m(2);
m(3), 0, -m(1);
-m(2), m(1), 0 ];
x=rand(3,4728);
B=A*x;
shows that it recovers the true, original m:
z = zeros(size(x,2),1);
Xbig = [ z x(3,:)' -x(2,:)';
-x(3,:)' z x(1,:)';
x(2,:)' -x(1,:)' z ];
Bbig = reshape(B',[],1);
m = Xbig\Bbig;
mtrue,m
mtrue = 3×1
0.3542 0.9438 0.0832
m = 3×1
0.3542 0.9438 0.0832

Bruno Luong on 7 Apr 2021
Edited: Bruno Luong on 7 Apr 2021
There might be a better mehod, at least more geometric, than linear system solving.
I understand you want to find m (3 x 1) such that
cross(m, x) = B.
m must be perpendicular to all columns of B, and B must be on a 2D plane.
If you make an SVD of B the smallest singular vector is then // to m, the two others form a basis of the plane where B live in.
You can project x and B on this plane, call the projection xp and Bp.
If you take their cross product
cross(xp,Bp)
you must find it equal to m*|xp|^2 (from triplet cross product properties).
So we can to estimate m simply as
m = mean( cross(xp,Bp) / |xp|^2).
There might be something similar with quaternion, all those are related somehow.
Test code
% Create fake test data
m = rand(3,1)
x = randn(3,10); % 10 must be replaced by 4728 in your case
B = cross(repmat(m,1,size(x,2)),x)
% Compute m from x and B
[U,~,~] = svd(B,0);
Q = U(:,1:2);
P = Q*Q';
Bp = P*B;
xp = P*x;
mrecover = mean(cross(xp,Bp,1)./sum(xp.^2,1),2)
Bruno Luong on 7 Apr 2021
Seem like a nice package Matt.

Matt J on 7 Apr 2021
Edited: Matt J on 7 Apr 2021
N=size(x,2);
C=func2mat( @(m)cross(repelem(m,1,N),x) , ones(3,1) , 'doSparse',0);
m=C\B(:);
A = [ 0, -m(3), m(2);
m(3), 0, -m(1);
-m(2), m(1), 0 ];

Jan Buchali on 8 Apr 2021
Unfortunaly nothing really worked out in my case. I also tried the lsqnonlin Solver but there I have other problems.
dipol = @(m,i) cross(m, x(i,:)') - B(:,i)';
m0 = [0.02; -0.01; 0];
lb = [-0.05; -0.05; -0.05];
ub = [0.05; 0.05; 0.05];
options.StepTolerance = 1e-10;
options.FunctionTolerance = 1e-10;
options.OptimalityTolerance = 1e-100;
options.MaxIterations = 1e3;
[m,resnorm,residual,exitflag,output] = lsqnonlin(@(m) dipol(m,i),m0,lb,ub,options);
In the upper case it says that: "the initial point is a local minimum", so that m is always m0. Can I do something on my options to optimize the solver?
Matt J on 8 Apr 2021
Also, if the columns of B are significantly non-coplanar, the solution will always be m=[0;0;0].