Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Matrix inverse in mathematica

If $A$ is an invertible $n \times n$ matrix, then $A\cdot A^{-1} = I$.

To get this statement in Mathematica, you need the assumption
MatrixPower[A, 0] = IdentityMatrix[n]

$Assumptions = {Element[A, Matrices[{n, n}]], 
                Det[A] != 0, ForAll[{A}, MatrixPower[A, 0] == IdentityMatrix[n]]}
TensorExpand[Inverse[A].A] // Simplify

Out[1]= IdentityMatrix[n].

Matrix as a tensor

A matrix can also be seen as a tensor of rank 2; i.e. a list with two levels, on for the columns and one for the rows. The tensor product $A \otimes A^{-1}$ corresponds in Mathematica with the outer product of the two lists. In this case it is a list with four levels, a four-tensor.

When matrices are viewed as tensors, the dot product is the same as a tensor product followed by a contraction. In Mathematica the matrix product $A\cdot A^{-1}$ can also be written as

TensorContract[Inverse[A]\[TensorProduct]A, {2, 3}] 

Question 1

How can I let Mathematica evaluate the above expression to IdentityMatrix[n]?

Involving more tensors

Suppose there is a second $n \times n$ matrix $B$. In that case one can think of more complex tensor products, for example $B \otimes A \otimes A^{-a}$. This is a rank 6 tensor.

Contracting slots 3 and 4 of the 6-tensor gives 4-tensor $B \otimes I$, where $I$ is the identity matrix. In Mathematica, this contraction can be written as

TensorContract[B\[TensorProduct]A\[TensorProduct]Inverse[A], {3, 4}] 

Question 2
Also for the above contraction, I would like Mathematica to use the idenity $A\cdot A^{-1} = I$. It should evaluate to

B\[TensorProduct]A\[TensorProduct]Inverse[A], {3, 4}] 

More over, I want Mathematica to use the same identity in more complicated tensor products, like $ B \otimes A \otimes B \otimes A^{-1} $. The last example should evaluate to

TensorTranspose[B \otimes I \otimes B, {{4, 5}]

How can this be realized with Mathematica 9?

Generalizations

Ihe above, there was allays only one invertible matrix which was called $A$.
Question 3
Can you take identity relations are into account for any invertible matrix, irrespectively of the names or number of invertible matrices?

share|improve this question
    
Maybe A\[TensorProduct]B ^:= IdentityMatrix[n] /; B == Inverse[A]? –  swish Feb 28 '13 at 14:33
    
@swish this may work in some cases, but not for all. Think of $A \otimes B \otimes A^{-1}$ for example. –  sjdh Feb 28 '13 at 16:18

1 Answer 1

up vote 5 down vote accepted
+50

Here's a rule that I think captures the spirit of what you're trying to do. (EDIT: Needed to shuffle some of the indices around to get the identity matrix into the right spot. SECOND EDIT: Treat the case of both left and right contraction.)

$Assumptions = 
 A ∈ Matrices[{n, n}] && 
  Inverse[A] ∈ Matrices[{n, n}] && Det[A] != 0 && 
  ForAll[{A}, MatrixPower[A, 0] == IdentityMatrix[n]] && 
  a ∈ Matrices[{n, n}] && b ∈ Matrices[{n, n}] && 
  c ∈ Matrices[{n, n}];

contractinv = 
 Quiet[{TensorContract[
     a___\[TensorProduct]A_\[TensorProduct]b___\[TensorProduct]B_\
\[TensorProduct]c___, {i1___, {m_, n_}, 
      i2___}] :> (TensorContract[
       TensorTranspose[
        a\[TensorProduct]IdentityMatrix[
          First[TensorDimensions[
            A]]]\[TensorProduct]b\[TensorProduct]c, 
        Join[Range[Min[m, n] - 1], {Max[m, n] - 1}, 
         Range[Min[m, n], Max[m, n] - 2]]], {i1, 
         i2} /. {i_ /; i > Max[m, n] :> i - 2}]) /; (B == Inverse[A] ||
         A == Inverse[B]) && m == 2 Length[{a}] + 2 && 
      n == 2 Length[{a, b}] + 3, 
   TensorContract[
     a___\[TensorProduct]A_\[TensorProduct]b___\[TensorProduct]B_\
\[TensorProduct]c___, {i1___, {m_, n_}, 
      i2___}] :> (TensorContract[
       TensorTranspose[
        a\[TensorProduct]IdentityMatrix[
          First[TensorDimensions[
            A]]]\[TensorProduct]b\[TensorProduct]c, 
        Join[Range[Min[m, n]], {Max[m, n] - 2}, 
         Range[Min[m, n] + 1, Max[m, n] - 3]]], {i1, 
         i2} /. {i_ /; i > Max[m, n] :> i - 2}]) /; (B == Inverse[A] ||
         A == Inverse[B]) && m == 2 Length[{a}] + 1 && 
      n == 2 Length[{a, b}] + 4}]

Which gives for a simple test:

TensorContract[
  a\[TensorProduct]Inverse[
    A]\[TensorProduct]b\[TensorProduct]A\[TensorProduct]c, {{2, 
    3}, {4, 7}, {8, 9}}] /. contractinv
(* TensorContract[a\[TensorProduct]b\[TensorProduct]c, {{2, 3}, {4, 7}}] *)

TensorContract[Inverse[A]\[TensorProduct]A,{{2,3}}]/.contractinv
(* IdentityMatrix[n] *)
share|improve this answer
    
This captures the spirit indeed. Can you explain how it works roughly? –  sjdh Feb 28 '13 at 19:32
    
This TensorContract[a\[TensorProduct]Inverse[A]\[TensorProduct]A, {{4, 5}}] /. contractinv gives $a$. It should give $a \otimes I$. Can you change that? –  sjdh Feb 28 '13 at 19:35
    
Roughly, it scans a tensor product inside a contraction to see if there are a matching pair of inverses that are being contracted. If so, it removes the pair, replaces them with an identity matrix, reshuffles the indices appropriately and then updates any remaining contraction indices. –  Xerxes Feb 28 '13 at 20:48
    
Nice solution @Xerxes. One detail, this example: TensorContract[Inverse[A]\[TensorProduct]B\[TensorProduct]b\[TensorProduct]A\[T‌​ensorProduct]c, {{2,7}, {4,10}, {8,9}}]/.contractinv gives TensorContract[Id\[TensorProduct]B\[TensorProduct]b\[TensorProduct]c, {{2,7}, {5,8}}] but it should be TensorContract[Id\[TensorProduct]B\[TensorProduct]b\[TensorProduct]c, {{2,7}, {4,8}}]. Any chance of correcting? Note I had to replace IdentityMatrix[3] with the symbol Id for version 9. Thanks. @sjdh –  Artur Gower Oct 30 at 13:52

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.