[Insight-users] Calculating least squares straight-line

[Sajendra]

Hi Matthias,

I'm not familiar with the itk FEM classes nor using them for data
fitting; however, If you just need simple linear least solution to an
overdetermined system can you not use the closed form solution?:
http://en.wikipedia.org/wiki/Linear_least_squares#Computation
http://mathworld.wolfram.com/LeastSquaresFitting.html

It is quite simple to compute. I recently hacked an ITK image filter
to perform that computation on input ITK image intensities so I could
use ITKs multithreading and IO. Let me know if that would help with
your problem.

Sajendra



On Tue, Nov 24, 2009 at 4:10 AM, Matthias Dodt
<matthias.dodt at mdc-berlin.de> wrote:
> Thanks for the hint - but i think the problem is that the matrix A has to be
> a square matrix. My system of equations is overdetermined, so maybe i should
> use a different solver...?
>
> I would need a MultipleValuedLinearOptimizer- but i am not sure which class
> to use-
>
> best,
>
> mat
>
> C.Cagatay Bilgin schrieb:
>>
>> I did not know this functionality existed in ITK and implemented
>> my own. Can you try running this on a smaller set of points and
>> provide those points to us so that I can give it a try as well.
>>
>> Also, I see the following lines in the source code that looks related
>> to your problem. If you have zero elements, just setting them to
>> something really small like 1.0e-16 might help, maybe.
>>
>> /* *******************************************************************
>> 00541    * FIX ME: itpack does not allow for any non-zero diagonal
>> elements
>> 00542    * so "very small" numbers are inserted to allow for a solution
>> 00543    *
>> 00544   int i;
>> 00545   doublereal fakeZero = 1.0e-16;
>> 00546
>> 00547   //insert "fake" zeros
>> 00548   for (i=0; i<static_cast<int>(m_Order); i++)
>> 00549   {
>> 00550     if ( (*m_Matrices)[0].Get(i,i) == 0.0)
>> 00551     {
>> 00552       (*m_Matrices)[0].Set(i,i,fakeZero);
>> 00553     }
>> 00554   }
>> 00555   // END FIXME
>> 00556
>>  *********************************************************************/
>>
>>
>> On Nov 23, 2009, at 10:22 AM, Matthias Dodt wrote:
>>
>>> Hi!
>>>
>>> I have some problems using the LinearSystemWrapperItpack class. I want to
>>> compute the straight-line wich fits a set of points best. Till now i got:
>>>
>>> void CalculateLeastSquares(){
>>>      using namespace fem;
>>>      //Least squares: b + A * x = y
>>>      //given Pi(xi,yi): 1b + xi * x = yi
>>>
>>>      fem::LinearSystemWrapperItpack lsw;
>>>      lsw.SetNumberOfMatrices(1u);
>>>      lsw.SetNumberOfVectors(1u);
>>>      lsw.SetNumberOfSolutions(1u);
>>>      lsw.SetSystemOrder(areaPoints->size());
>>>      lsw.SetMaximumNonZeroValuesInMatrix(areaPoints->size() *4);
>>>      lsw.InitializeMatrix(0);
>>>      lsw.InitializeSolution(0);
>>>      lsw.InitializeVector(0);
>>>
>>>
>>>      //Matrix<double,2> coeff;
>>>      //Vector<double,1> y;
>>>      for(unsigned int t=0;t<areaPoints->size();++t){
>>>          //coeff(0,i) = 1;
>>>          std::cout << "Setting: " << t << "," << areaPoints->at(t)[0] <<
>>> "," << areaPoints->at(t)[1] << std::endl;
>>>          lsw.SetMatrixValue(t,0u,1u,0);
>>>
>>>  lsw.SetMatrixValue(t,1,static_cast<float>(areaPoints->at(t)[0]),0);
>>>          //coeff(1,i) = static_cast<double>(areaPoints->at(i)[0]);
>>>
>>>  lsw.SetSolutionValue(t,static_cast<float>(areaPoints->at(t)[1]),0);
>>>          //y(i) = static_cast<double>(areaPoints->at(i)[1]);
>>>      }
>>>
>>>      lsw.Solve();
>>>
>>>      for(unsigned int i=0;i<lsw.GetNumberOfSolutions();++i){
>>>          std::cout << "Coefficient: " << lsw.GetSolutionValue(i,0);
>>>      }
>>>
>>>
>>>
>>>  }
>>>
>>> Anyway the .Solve() crashes with:
>>>
>>> itk::FEMExceptionItpackSolver (0x66c1a0)
>>> Location: "LinearSystemWrapperItpack::Solve"
>>> File:
>>> /Users/mat/Documents/InsightToolkit-3.14.0/Code/Numerics/FEM/itkFEMLinearSystemWrapperItpack.cxx
>>> Line: 659
>>> Description: Error: A diagonal element is not positive
>>>
>>> Somebody has a clue?
>>>
>>> thanks!
>>>
>>> _____________________________________
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at
>>> http://www.kitware.com/opensource/opensource.html
>>>
>>> Kitware offers ITK Training Courses, for more information visit:
>>> http://www.kitware.com/products/protraining.html
>>>
>>> Please keep messages on-topic and check the ITK FAQ at:
>>> http://www.itk.org/Wiki/ITK_FAQ
>>>
>>> Follow this link to subscribe/unsubscribe:
>>> http://www.itk.org/mailman/listinfo/insight-users
>>
>
> --
> ------------------------------------------------
> Matthias Dodt
>
> Scientific Programmer at Bioinformaitcs platform AG Dieterich
>
> Berlin Institute for Medical Systems Biology at the
> Max-Delbrueck-Center for Molecular Medicine
> Robert-Roessle-Strasse 10, 13125 Berlin, Germany
>
> fon: +49 30 9406 4261
> email: matthias.dodt at mdc-berlin.de
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>