[Insight-users] Calculating least squares straight-line

[Matthias Dodt]

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