opto abhi
opto abhi

Reputation: 317

matrix multiplication error between MATLAB and .NET

Note: this question is related to my previous question here on Matrix right division

When I compare the final result of the following example in both MATLAB and C#, I notice that there is a noticeable difference. Why is this so ?

The result of finding matrix inverses seems to tally, but the A*A^-1 seems to be way off.

Example in MATLAB:

>> a = [1 2 3; 4 5 6; 7 8 10]

a =

     1     2     3
     4     5     6
     7     8    10
>> inv(a)

ans =

   -0.6667   -1.3333    1.0000
   -0.6667    3.6667   -2.0000
    1.0000   -2.0000    1.0000
>> z = mtimes(a, inv(a))
>> z

z =

   1.0000e+00  -4.4409e-16  -1.1102e-16
   1.3323e-15   1.0000e+00  -2.2204e-16
   2.2204e-15  -2.6645e-15   1.0000e+00

Same data in C#: using CSML Matrix Library

//using CSML Matrix Library
public static Matrix operator *(Matrix A, Matrix B)
{
    if (A.ColumnCount != B.RowCount)
    throw new ArgumentException("Inner matrix dimensions must agree.");
    Matrix C = new Matrix(A.RowCount, B.ColumnCount);

        for (int i = 1; i <= A.RowCount; i++)
        {
            for (int j = 1; j <= B.ColumnCount; j++)
            {
                C[i, j] = Dot(A.Row(i), B.Column(j));
            }
        }

    return C;
}
Console.WriteLine(e1);
1;      2;      3;      \
4;      5;      6;      \
7;      8;      10;     \

Console.WriteLine(e1.Inverse());
-0.666666666666667;     -1.33333333333333;      1;      \
-0.666666666666669;     3.66666666666667;       -2;     \
1;      -2;     1;      \

Console.WriteLine(e1 * e1.Inverse());
0.999999999999999;      1.77635683940025E-15;   -8.88178419700125E-16;  \
-5.32907051820075E-15;  1.00000000000001;       -3.5527136788005E-15;   \
-1.06581410364015E-14;  3.5527136788005E-15;    0.999999999999996;      \

Upvotes: 1

Views: 128

Answers (1)

Matthew Thirkettle
Matthew Thirkettle

Reputation: 442

Both results seem sensible. MATLAB computes inv(A) using row-reduction. Each step of row-reduction results in numerical error (MATLAB interprets 1/3 as a decimal with a finite number of decimal places). So, due to numerical error, I would expect the elements of inv(A) to be off by 10^{-16}. The result that

a*inv(a) = 1 +/- 10^{-16} along the diagonals and

a*inv(a) = +/- 10^{-16} along the off-diagonals

is consistent with a*inv(a) equalling the identity matrix plus or minus some numerical error.

Upvotes: 2

Related Questions