Ask Your Question
1

Bug in matrix multiplication and inversion?

asked 2016-04-18 16:06:07 +0200

Thomas gravatar image

I'm using sage 7.1 and this code here returns false.

M = matrix ([[sqrt(1-1/2),1/2],[1,-sqrt(1-1/2)]])
N = matrix ([[sqrt(1-1/3),1],[1/3,-sqrt(1-1/3)]])
(M*N).inverse() == N.inverse()*M.inverse()

What is going on?

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
1

answered 2016-04-18 17:16:10 +0200

tmonteil gravatar image

If E denotes the difference;

sage: E = (M*N).inverse() - N.inverse()*M.inverse()
sage: E
[-sqrt(2/3)*sqrt(1/2) + 6/(6*sqrt(2/3)*sqrt(1/2) + 1) - 6*(3*sqrt(2/3) - sqrt(1/2))*(sqrt(2/3) - 2*sqrt(1/2))/((6*sqrt(2/3)*sqrt(1/2) + 1)^2*(sqrt(2/3)*sqrt(1/2) + (3*sqrt(2/3) - sqrt(1/2))*(sqrt(2/3) - 2*sqrt(1/2))/(6*sqrt(2/3)*sqrt(1/2) + 1) + 1)) - 1                                                           -1/2*sqrt(2/3) + sqrt(1/2) + 3*(sqrt(2/3) - 2*sqrt(1/2))/((6*sqrt(2/3)*sqrt(1/2) + 1)*(sqrt(2/3)*sqrt(1/2) + (3*sqrt(2/3) - sqrt(1/2))*(sqrt(2/3) - 2*sqrt(1/2))/(6*sqrt(2/3)*sqrt(1/2) + 1) + 1))]
[                                                           sqrt(2/3) - 1/3*sqrt(1/2) - 2*(3*sqrt(2/3) - sqrt(1/2))/((6*sqrt(2/3)*sqrt(1/2) + 1)*(sqrt(2/3)*sqrt(1/2) + (3*sqrt(2/3) - sqrt(1/2))*(sqrt(2/3) - 2*sqrt(1/2))/(6*sqrt(2/3)*sqrt(1/2) + 1) + 1))                                                                                                                   -sqrt(2/3)*sqrt(1/2) + 1/(sqrt(2/3)*sqrt(1/2) + (3*sqrt(2/3) - sqrt(1/2))*(sqrt(2/3) - 2*sqrt(1/2))/(6*sqrt(2/3)*sqrt(1/2) + 1) + 1) - 1/6]

You get a matrix whose entries are symbolic expressions. If we write:

sage: E == 0
False

We get false because the entries of E are not all zero, as symbolic expressions, but they are actually zero if we simplify all entries one by one:

sage: E[0][0]
-sqrt(2/3)*sqrt(1/2) + 6/(6*sqrt(2/3)*sqrt(1/2) + 1) - 6*(3*sqrt(2/3) - sqrt(1/2))*(sqrt(2/3) - 2*sqrt(1/2))/((6*sqrt(2/3)*sqrt(1/2) + 1)^2*(sqrt(2/3)*sqrt(1/2) + (3*sqrt(2/3) - sqrt(1/2))*(sqrt(2/3) - 2*sqrt(1/2))/(6*sqrt(2/3)*sqrt(1/2) + 1) + 1)) - 1
sage: E[0][0].full_simplify()
0
sage: E[0][1].full_simplify()
0
sage: E[1][1].full_simplify()
0
sage: E[1][0].full_simplify()
0

Note that checking whether a symbolic expression is equal to zero is undecidable in general (though we should admit that in the present case, this is just because symbolics is weak in Sage). Since your symbolic expressions represent algebraic numbers, you can change the ring over which E is defined to be the field of algebraic numbers. In this field, the equalty is decidable, and Sage knows how to handle it well:

sage: E.change_ring(QQbar) == 0
True
edit flag offensive delete link more

Comments

2

Not that one can also simplify the matrix E directly:

sage: E.simplify_full()
[0 0]
[0 0]

In terms of weakness of symbolics in Sage, is there a reason for not trying to apply simplify_full when one tests an expression for zero? (Note that I am very much illiterate regarding Sage's symbolics!)

B r u n o gravatar imageB r u n o ( 2016-04-18 18:17:27 +0200 )edit

Damn, i was looking for full_simplify() but tab completion did not gave anything !

tmonteil gravatar imagetmonteil ( 2016-04-18 21:41:53 +0200 )edit

OK thanks. I did actually try simplify on the difference, but not simplify_full()

Thomas gravatar imageThomas ( 2016-04-19 12:57:15 +0200 )edit

... and by the way, Mathematica just does it.

Thomas gravatar imageThomas ( 2016-04-19 12:57:51 +0200 )edit
2

As a general rule, when you are doing purely algebraic computation, an advice is to tell Sage that you're doing algebra, and not deal with symbolics. Here, you can do:

sage: M = matrix (QQbar,[[sqrt(1-1/2),1/2],[1,-sqrt(1-1/2)]])
sage: N = matrix (QQbar,[[sqrt(1-1/3),1],[1/3,-sqrt(1-1/3)]])
sage: (M*N).inverse() == N.inverse()*M.inverse()
True

And another solution is not to work on QQbar but staying in the symbolic world and using the much more powerful method is_zero:

sage: M = matrix ([[sqrt(1-1/2),1/2],[1,-sqrt(1-1/2)]])
sage: N = matrix ([[sqrt(1-1/3),1],[1/3,-sqrt(1-1/3)]])
sage: ((M*N).inverse() - N.inverse()*M.inverse()).is_zero()
True
B r u n o gravatar imageB r u n o ( 2016-04-19 14:55:52 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2016-04-18 16:06:07 +0200

Seen: 414 times

Last updated: Apr 18 '16