forum.alglib.net

ALGLIB forum
It is currently Tue Dec 10, 2024 11:02 am

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 1 post ] 
Author Message
 Post subject: Matrix inversion produces unexpected results
PostPosted: Sun Jun 05, 2022 9:06 pm 
Offline

Joined: Sun Jun 05, 2022 8:36 pm
Posts: 1
Please see the first code excerpt below along with an input matrix, Z, and an output matrix, Zinverse. Zinverse is produced by the alglib.rmatrixinverse; and it is obviously wrong. However, in the second code excerpt, alglib.rmatrixinverse produces a correct result. I do not see reason for the difference in the results. As far as I can tell, the problem is occuring within alglib.rmatrixinverse. Can you please help me understand (1) how this results difference is happening, and (2) how to diagnose issues within alglib methods?

Thanks in advance.





// search for the non-singular scaling matrix with largest termCount
termCount = 6; // (m / 4) + 2;
while (searching && termCount >= MIN_TERM_COUNT)
{
matrix.Multiply(Yt, Y, termCount, ref Z);
if (alglib.rmatrixdet(Z, termCount) != 0)
{
File.WriteAllText(Path.Combine(dataPath, $"Z_{termCount}.txt"), matrix.FormatMatrix(Z, "Z", fieldWidth));
searching = false;
// compute inverse
alglib.rmatrixinverse(ref Z, termCount, out aglibInfo, out aglibReport); // inverse returns in Z matrix
if (aglibInfo == 1) // inversion was successful
{
File.WriteAllText(Path.Combine(dataPath, $"Zinverse_{termCount}.txt"), matrix.FormatMatrix(Z, "Zinverse", fieldWidth));

matrix.Multiply(Z, Yt, ref S);
matrix.Multiply(S, x, ref a);
SetScalingConstants(a);
TermCount = termCount;
rv = true;
break;
}
else
throw new InvalidOperationException($"alglib.rmatrixinverse: aglibInfo = {aglibInfo}");
}
else
{
matrix.Initialize(ref Z, Matrix.InitializationTypes.Null);
termCount--;
}
}

File.WriteAllText(Path.Combine(dataPath, $"Y_{termCount}.txt"), matrix.FormatMatrix(Y, "Y", fieldWidth));
File.WriteAllText(Path.Combine(dataPath, $"Yt_{termCount}.txt"), matrix.FormatMatrix(Yt, "Yt", fieldWidth));
File.WriteAllText(Path.Combine(dataPath, $"S_{termCount}.txt"), matrix.FormatMatrix(S, "S", fieldWidth));


Z = { 6.00 -9.28 3.21 -1.81 -0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ -9.28 17.65 -6.75 3.21 0.01 -0.04 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 3.21 -6.75 2.69 -1.19 -0.01 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ -1.81 3.21 -1.19 0.60 0.00 -0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ -0.01 0.01 -0.01 0.00 0.00 -0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.01 -0.04 0.02 -0.01 -0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }

Zinverse = { 10646.89 201458.61 164016.83 -660336.75 -4243385.60 4725122.52 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 201458.61 3918229.40 3163567.51 -12893771.96 -89957384.22 90987096.05 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 164016.83 3163567.51 2560885.83 -10397909.15 -70690760.75 73714030.96 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ -660336.75 -12893771.96 -10397909.15 42453443.94 299596879.81 -298962625.64 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ -4243385.60 -89957384.22 -70690760.75 299596879.82 2674360516.42 -2003109303.81 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 4725122.52 90987096.05 73714030.96 -298962625.64 -2003109303.73 2126085570.33 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }
{ 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 }


public void Test_ZInverse()
{
int info;
alglib.matinvreport rep;

Matrix matrix = new Matrix();
double[,] I;
double[,] P = null;
int termCount = 6;
double[,] X= null;
double[,] Z = null;
double[,] Zinverse = null;

initializeZInputSet(ref Zinverse, ref X);
if (alglib.rmatrixdet(Zinverse, termCount) != 0)
{
//lglib.rmatrixinverse(ref Zinverse, termCount, out info, out rep);
alglib.rmatrixinverse(ref Zinverse, termCount, out info, out rep);
Assert.True(info == 1);
}
else
Assert.True(false);

initializeZInputSet(ref Z, ref X);
matrix.Multiply(Z, Zinverse, termCount, ref P);
I = matrix.NewMatrix(termCount, termCount, Matrix.InitializationTypes.Identity);
Assert.True(matrix.IsEqual(P, I));
}


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 1 post ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 2 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group