Message Boards Message Boards

0
|
6762 Views
|
8 Replies
|
4 Total Likes
View groups...
Share
Share this post:

Numerical Error with Matrix operations

Posted 10 years ago

A is a 3x3 matrix, b is a 3x1 vector. I try to convert [A|b], a 3x4 matrix, to [I|0]. So the formula is right multiple [A|b] by a matrix:

$$ \begin{bmatrix} A^{-1} & -A^{-1}b\ 0 & 1\ \end{bmatrix} $$

Here's the code in Mathematica 10

A=Table[RandomReal[],{i,3},{j,3}]

b=Table[RandomReal[],{i,3}]

M=Transpose[Join[A, {b}]]

ret=M.Join[Transpose[Join[Inverse[A], {-Inverse[A].b}]], {{0, 0, 0, 1}}]

The last column of ret is always non-zero. The ret will often look like for example

$$ \begin{bmatrix} 1. & 0. & 0. & -0.979237\ 0. & 1. & 4.44089*10^{-16} & 2.17005\ -8.8817810^{-16} & -5.5511210^{-17} & 1. & -1.79119 \end{bmatrix} $$

How can I make Mathematica do it right?

POSTED BY: Gordon
8 Replies

A straightfoward method to check the matrix product you asked for: define your matrix by blocks using ArrayFlatten:

A = Table[RandomReal[], {i, 3}, {j, 3}];
b = List /@ Table[RandomReal[], {i, 3}]; (*So it is a matrix*)
M = ArrayFlatten[{{A, b}}];
ret = ArrayFlatten[{{Inverse[A], -Inverse[A].b}, {0, 1}}]; (*Notice numbers represents matrices for ArrayFlatten*)
M.ret // Chop // MatrixForm

Which returns the 3x4 identity matrix we were expecting.

Posted 10 years ago

It's so good to learn ArrayFlatten. Helps a lot.

Thank you

POSTED BY: Gordon

I guess I should state that Mathematica is doing it right. Truncation error is not avoidable with finite precision arithmetic. As for the problem at hand, there are a couple of things that seem to be amiss.

(1) Create the augmented matrix {A|b] correctly. This can be done with M = Join[A, List /@ b, 2]. Your version transposes A in the process.

(2) Now just use RowReduce on it. No inverses, no multiplications.

Here is an example.

SeedRandom[1111];
aA = Table[RandomReal[], {i, 3}, {j, 3}]
b = Table[RandomReal[], {i, 3}]
mat = Join[aA, List /@ b, 2]
rred = RowReduce[mat]

Out[434]= {{0.0777960697813, 0.556265936296, 
  0.807681686834}, {0.527982841158, 0.555876762387, 
  0.981090055861}, {0.907596678267, 0.364390115148, 0.722450567448}}

Out[435]= {0.420186553434, 0.1560084863, 0.13002033209}

Out[436]= {{0.0777960697813, 0.556265936296, 0.807681686834, 
  0.420186553434}, {0.527982841158, 0.555876762387, 0.981090055861, 
  0.1560084863}, {0.907596678267, 0.364390115148, 0.722450567448, 
  0.13002033209}}

Out[437]= {{1, 0., 0., 0.415692234671}, {0, 1, 0., 4.46163163168}, {0, 0, 1, -2.59261340179}}

Check that this agrees with LinearSolve (I assume that was the purpose).

rred[[All, -1]] == LinearSolve[aA, b]

Out[438]= True
POSTED BY: Daniel Lichtblau
Posted 10 years ago

Thanks for pointing out my mistake on [A|b] construction. That's the best answer.

POSTED BY: Gordon

My apologies, I forgot about the last column as the title mentioned 'numerical errors'. Thanks for completing the answer.

Gordon,

Replace -Inverse[A].b with -b.Inverse[A] and then use Chop as suggested by Guillermo.

A = Table[RandomReal[], {i, 3}, {j, 3}] ;
b = Table[RandomReal[], {i, 3}] ;
M = Transpose[Join[A, {b}]] ;
ret = M.Join[
    Transpose[Join[Inverse[A], {-b.Inverse[A]}]], {{0, 0, 0, 1}}] ;
ret // MatrixForm
ret // Chop // MatrixForm

I.M.

POSTED BY: Ivan Morozov
Posted 10 years ago

Wow, it works. Why? What's the magic behind that?

POSTED BY: Gordon

Remember that when you perform numerical computations with floating point numbers they have a limited precision, which might lead to numerical errors. It is usual to assume in such cases that results that are very close to zero are zero indeed. Mathematica's Chop function can be used for that. Just try:

Chop[ret]

You can also check it is a precision error by using arbitray precision for your numers:

A=Table[RandomReal[WorkingPrecission->30],{i,3},{j,3}]
b=Table[RandomReal[WorkingPrecission->30],{i,3}]
M=Transpose[Join[A, {b}]]
ret=M.Join[Transpose[Join[Inverse[A], {-Inverse[A].b}]], {{0, 0, 0, 1}}]

Which will make those almost-zeros even smaller.

You can also check this tutorial for more information.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract