1

I have a linear algebraic equation of the form Ax=By. Where A is a matrix of 6x5, x is vector of size 5, B a matrix of 6x6 and y vector of size 6. A, B and y are known variables and their values are accessed in real time coming from the sensors. x is unknown and has to find. One solution is to find Least Square Estimation that is x = [(A^T*A)^-1]*(A^T)B*y. This is conventional solution of linear algebraic equations. I used Eigen QR Decomposition to solve this as below

matrixA = getMatrixA();
matrixB = getMatrixB();
vectorY = getVectorY();
//LSE Solution
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> dec1(matrixA);
vectorX = dec1.solve(matrixB*vectorY);//

Everything is fine until now. But when I check the errore = Ax-By, its not zero always. Error is not very big but even not ignorable. Is there any other type of decomposition which is more reliable? I have gone through one of the page but could not understand the meaning or how to implement this. Below are lines from the reference how to solve the problem. Could anybody suggest me how to implement this?

The solution of such equations Ax = Byis obtained by forming the error vector e = Ax-By and the finding the unknown vector x that minimizes the weighted error (e^T*W*e), where W is a weighting matrix. For simplicity, this weighting matrix is chosen to be of the form W = K*S, where S is a constant diagonal scaling matrix, and K is scalar weight. Hence the solution to the equation becomes 

x = [(A^T*W*A)^-1]*(A^T)*W*B*y

I did not understand how to form the matrix W.

Wafeeq
  • 869
  • 1
  • 16
  • 34
  • You can try `Eigen::FullPivLU`, it might be a tiny bit more accurate. It will not solve your underlying cause though: You're solving the normal equations, hence a least squares approximation of your original system. – Ela782 Jan 27 '15 at 11:45

1 Answers1

2

Your statement " But when I check the error e = Ax-By, its not zero always. " almost always will be true, regardless of your technique, or what weighting you choose. When you have an over-described system, you are basically trying to fit a straight line to a slew of points. Unless, by chance, all the points can be placed exactly on a single perfectly straight line, there will be some error. So no matter what technique you use to choose the line, (weights and so on) you will always have some error if the points are not colinear. The alternative would be to use some kind of spline, or in higher dimensions to allow for warping. In those cases, you can choose to fit all the points exactly to a more complicated shape, and hence result with 0 error.

enter image description here

So the choice of a weight matrix simply changes which straight line you will use by giving each point a slightly different weight. So it will not ever completely remove the error. But if you had a few particular points that you care more about than the others, you can give the error on those points higher weight when choosing the least square error fit.

For spline fitting see:

http://en.wikipedia.org/wiki/Spline_interpolation

For the really nicest spline curve interpolation you can use Centripital Catmull-Rom, which in addition to finding a curve to fit all the points, will prevent unnecessary loops and self intersections that can sometimes come up during abrupt changes in the data direction.

Catmull-rom curve with no cusps and no self-intersections

Community
  • 1
  • 1
Ted
  • 3,212
  • 25
  • 20
  • It's called "Least Squares Approximation" for a reason: it is an approximation. You are trying to find a line of best fit. If you want zero error, you could always just using cubic splines to fit every data point, but now you have a complicated polynomial rather than a straight line, and it doesn't necessarily fix anything for you. – Cloud Mar 28 '14 at 17:14
  • @Ted So you mean I need a curve fitting technique to minimize the error? Normally what I know about curve fitting techniques is that we have some random data not the linear one a equation of some degree can be extracted using some Matlab tool etc I remember from my past experience I did something like this. But today in this example I don't have a past data.I have only a one instant point of value. How could I know if it fits or not? – Wafeeq Mar 28 '14 at 17:28
  • @Dogbert how to implement `cubic spline` ? What do you mean I did not understand you. – Wafeeq Mar 28 '14 at 17:30
  • I don't understand how you can have both an over-described system and only one point? – Ted Mar 28 '14 at 17:30
  • Use the same matrix values you were using for your A and B, but just fit them to a spline. – Ted Mar 28 '14 at 17:31
  • @Ted At one instant I have output a vector `x` of `5` values as I have described in my question. – Wafeeq Mar 28 '14 at 17:34
  • Ax=By is the linear equation you built for the linear case. But with a curve, you will have higher order terms. This gets kind of confusing when working with Matrices and vectors, but the basic premise still holds. Your equation needs to probably have a polynomial with the same number of unknowns as points you are trying to fit. If that is five then it becomes: Y = Ax^4 + Bx^3 + Cx^2 + Dx + E. – Ted Mar 28 '14 at 17:43
  • You are no longer looking at linear algebra at all, but rather at finding the exact solution for the polynomial that fits all your points. – Ted Mar 28 '14 at 17:44
  • I have added some references to some spline fitting approaches. – Ted Mar 28 '14 at 17:53
  • @Ted I had a look into your references. Actually previously I had not any idea about splines. I have gone through several tutorials / videos to understand the basic idea of the splines. Let me explain you what did I understand. For example we have a series of data `[(1,1),(2,2),(3,3),....;(n,n)]` just a linear data. If the data is not linear or sometimes a noise comes out of the sensors then we can use spline interpolation to get a smoother one or just the output where we don't have any control point. Or I have a function `F(x) = x²`, I want that for each value of `x` , `F(x)-x² ==0` – Wafeeq Apr 02 '14 at 16:59
  • @Ted In my case. I have variables `A`,`B` and `y` which are continuously changing at new instant. For example at `t = 10`, my `A`,`B` and `y` are different and output x calculated by `QR decomposition` would be different from everything calculated at `t = 9`. I am using the `x` calculated at `t = 10` (current time). And if I check my solution `Ax-By` at `t=9` its not complete zero sometimes. and again if I check my solution at `t=10` gives some error. I don't need to find the values in between `t=9` and `t=10`. I just want that `Ax-By==0` always. `splines` methods are applied on four points. – Wafeeq Apr 02 '14 at 17:07