Calculating a matrix with a specific form

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]

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?"

Answer 1

[James Tursa]

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).

[Jan Buchali]

No, m is the magnetic dipol, I have to estimate. X is the magnetic field and B the residual torque in X,Y and Z axis. Because one column of the matrix is one data point i have to use A = B*X'*inv(X*X') to solve it, that´s why your code is not useable for me.

[Matt J]

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

Answer 2

[Bruno Luong]

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)

[Matt J]

This is essentially done for you in planarFit (Download),

fobj=planarFit(B);

m=fobj.normal(:);
c=cross(m,x(:,1));

m=m*norm(B(:,1))/norm(c)*sign(c.'*B(:,1));

[Bruno Luong]

Seem like a nice package Matt.

Answer 3

[Matt J]

Another solution, using func2mat (Download).

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  ];

Answer 4

[Jan Buchali]

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;

options = optimset('Diagnostics','off','Display','final','GradObj','on');   

[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]

An iterative solution is overkill for this problem, which has an analytical solution.

We need to understand why you are not accepting one of the other solutions that have already been given to you. They have all been tested and found to work.

[Bruno Luong]

"... 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."

"dipol = @(m,i) cross(m, x(i,:)') - B(:,i)';"

• What is "i" in you script?
• There must be something that is not coherent in the dimension of the matrices or the way you setup lsqnonlin problem.

This show all method give the same result:

% Create fake test data
m = rand(3,1)
p = 4728;
x = randn(3,p); % 10 must be replaced by 4728 in your case
B = cross(repmat(m,1,size(x,2)),x);

% Compute m from x and B using Bruno's method
[U,~,~] = svd(B,0);
Q = U(:,1:2);
P = Q*Q';
Bp = P*B;
xp = P*x;
mmBruno = mean(cross(xp,Bp,1)./sum(xp.^2,1),2)

% Compute m from x and B using James's method
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);
mJames = Xbig\Bbig

% Compute m from x and B using lsqnonlin
dipol = @(m) cross(repmat(m,1,size(x,2)),x) - B;
m0 = zeros(3,1);
options = optimset('Diagnostics','off','Display','final','GradObj','on');   
mlsqnonlin = lsqnonlin(dipol,m0,[],[],options)

[Matt J]

This shows that all 5 methods discussed produce the same result,

% Create fake test data
m = rand(3,1)

m = 3×1

0.3873 0.7200 0.3608

p = 4728;
x = randn(3,p); % 10 must be replaced by 4728 in your case
B = cross(repmat(m,1,size(x,2)),x);



% Compute m from x and B using Bruno's method
[U,~,~] = svd(B,0);
Q = U(:,1:2);
P = Q*Q';
Bp = P*B;
xp = P*x;
mmBruno = mean(cross(xp,Bp,1)./sum(xp.^2,1),2)

mmBruno = 3×1

0.3873 0.7200 0.3608

% Compute m from x and B using James's method
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);
mJames = Xbig\Bbig

mJames = 3×1

0.3873 0.7200 0.3608

% Compute m from x and B using FEX submission func2mat()
N=size(x,2);
C=func2mat( @(m)cross(repelem(m,1,N),x) , ones(3,1) , 'doSparse',0);
mMatt1=C\B(:)

mMatt1 = 3×1

0.3873 0.7200 0.3608

% Compute m from x and B using FEX submission planarFit()
fobj=planarFit(B);
m=fobj.normal(:);
c=cross(m,x(:,1));
mMatt2=m*norm(B(:,1))/norm(c)*sign(c.'*B(:,1))

mMatt2 = 3×1

0.3873 0.7200 0.3608

% Compute m from x and B using lsqnonlin
M = repmat(m,1,size(x,2));
dipol = @(m) cross(repmat(m,1,size(x,2)),x) - B;
m0 = zeros(3,1);
options = optimset('Diagnostics','off','Display','final','GradObj','on');   
mlsqnonlin = lsqnonlin(dipol,m0,[],[],options)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

mlsqnonlin = 3×1

0.3873 0.7200 0.3608

[James Tursa]

I wonder if there are observability issues with the actual data OP has, and maybe this is leading to full rank related problems downstream?

[Matt J]

Also, if the columns of B are significantly non-coplanar, the solution will always be m=[0;0;0].