Lsqnonlin problem for calibrating accelerometer

조회 수: 1 (최근 30일)
raphael ricard
raphael ricard 2021년 7월 1일
편집: John D'Errico 2021년 7월 1일
Hi, i need to calibrate an accelerometer with a method that i found in a paper. I linked a photo where you can see the method they propose. I have a bank of data that correspond to the matrix Xsense that is filled with values of the readings in [ x y z ]. In the paper they suggest using a non linear least square curve fitting to fit the values of the accelerometer on an ellipse. In the paper, they say the we have to find Q which is a positive definite matrix and b, a vector of the coordinates for the center [cx cy cz]. To find these values, they suggest using the lsqnonlin of matlab with the function (transpose(xsense − b) * Q *(xsense − b) - 1 ). From what I understood, the lsqnonlin needs to have a start point x0 and in all the example i have seen, x0 is a vector linked to the output of the fonction. My problem is that Q is a definite matrix that i dont know the dimensions and b is a vector 1x3. How can i use lsqnonlin to find Q and b ? I can't find a way to fit the initial condition of Q and b in the x0 structure needed by lsqnonlin.

답변 (1개)

John D'Errico
John D'Errico 2021년 7월 1일
편집: John D'Errico 2021년 7월 1일
Um, the solution is not to find Q. Because it is difficult to constrain a matrix to be positive definite. You don't want to try to do that using lsqnonlin (regardless of what they say.)
Instead, you want to find a triangular matrix, I'll call it R, such that R is UPPER triangular. Then compute Q from R.
Q = R' * R
That is, irregardless what the matrix R is, The product, R' * R will ALWAYS be at least positive semi-definite. The semi-definite case comes from when R is itself singular, and you can even prevent that from happening, if you are careful.
So, for example, pick ANY matrix R.
format long g
R = triu(randn(3))
R = 3×3
-1.73472568671219 0.133100695588286 0.421640106149951 0 2.10509564761744 -0.485099810835556 0 0 -1.18054675031493
Is Q positive semi-definite?
Q = R'*R
Q = 3×3
3.00927320813907 -0.230893195556259 -0.731429922686373 -0.230893195556259 4.449143480784 -0.965060909033498 -0.731429922686373 -0.965060909033498 1.80679283526597
It is by construction symmetric. The eigenvalues are...
eig(Q)
ans = 3×1
1.17173143780017 3.32936783906208 4.76411024732678
SURPRISE! It works! Actually, the classic test in MATLAB to see if a matrix is SPD, is to apply chol to it. In fact, what we did was we created Q from the cholesky factorization. So chol should simply recover R. Did it?
chol(Q)
ans = 3×3
1.73472568671219 -0.133100695588286 -0.421640106149951 0 2.10509564761744 -0.485099810835556 0 0 1.18054675031493
Well, it did, but for the signs on some of those elements. That is arbitrary. This is like a sqrt, the sign is indeterminate.
The point is, you need to search for only a upper triangular matrix R that does what you need. You could also solve this for a LOWER triangular matrix. But then, I suppose you need to call it L. :) Inside the objective function, form Q from R. This will be sufficient to span the set of positive definite matrices that you will need to search over. And best of all, this will be eminently solvable using lsqnonlin. You will search over UPPER triangular matrices R. So instead of 9 elements, you need only search over 6 unknowns.
When will there be a problem, and Q is only semi-definite? Again, that happens only when one of the diagonal elements of R is sufficiently small in magnitude compared to the other diagonal elements. You will need to watch for that, if you care that Q must always be strictly positive definite. The case where Q becomes singular from this is equivalent to an ellipsoid that contains zero volume. So effectively, you have an ellipsoid that has compressed into a plane, or worse, into a stright line. Those are the degenerate cases. But your data should make it painfully obvious when that would happen.

카테고리

Help CenterFile Exchange에서 Operating on Diagonal Matrices에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by