7
$\begingroup$

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?

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

1 Answer 1

7
+50
$\begingroup$

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] *)
$\endgroup$
4
  • $\begingroup$ This captures the spirit indeed. Can you explain how it works roughly? $\endgroup$
    – sjdh
    Feb 28, 2013 at 19:32
  • $\begingroup$ This TensorContract[a\[TensorProduct]Inverse[A]\[TensorProduct]A, {{4, 5}}] /. contractinv gives $a$. It should give $a \otimes I$. Can you change that? $\endgroup$
    – sjdh
    Feb 28, 2013 at 19:35
  • $\begingroup$ 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. $\endgroup$
    – Xerxes
    Feb 28, 2013 at 20:48
  • $\begingroup$ Nice solution @Xerxes. One detail, this example: TensorContract[Inverse[A]\[TensorProduct]B\[TensorProduct]b\[TensorProduct]A\[TensorProduct]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 $\endgroup$
    – Art Gower
    Oct 30, 2014 at 13:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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