<h1>Matrices</h1>

<p><a href="https://docs.sympy.org/latest/tutorial/matrices.html">From</a></p>


```julia
    >>> from sympy import *
    >>> init_printing(use_unicode=True)
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>In <code>SymPy</code>, matrices can be store using <code>Julia</code>'s <em>generic</em> <code>Matrix&#123;T&#125;</code> type where <code>T &lt;: Sym</code> <em>or</em> using SymPy's matrix type, wrapped in a <code>SymMatrix</code> type by <code>SymPy</code>. This tutorial shows how to use the underlying <code>SymMatrix</code> values. To construct a matrix of symbolic values is identical to construction a matrix of numeric values within <code>Julia</code>, and will be illustrated at the end.</p>
</li>
</ul>

<ul>
<li><p>methods for <code>SymMatrix</code> objects use the dot call syntax. As a convenience, this will also work for <code>Array&#123;Sym&#125;</code> objects. The returned value may be a <code>SymMatrix</code>, not an <code>Array&#123;Sym&#125;</code>.</p>
</li>
<li><p>The matrix constructor in SymPy using a vector of row vectors does <em>not</em> work in <code>SymPy</code>, as of newer versions (it does not work with version 1.5.1 of sympy and PyCall). This style is used in this document. The user of <code>SymPy</code> can eaesily avoid this specification, using Julia's matrix construction techniques.</p>
</li>
</ul>


```julia
using SymPy
```

<hr />

<p>To make a matrix in SymPy, use the <code>Matrix</code> object.  A matrix is constructed by providing a list of row vectors that make up the matrix.  For example, to construct the matrix</p>


$$
   \left[\begin{array}{cc}1 & -1\\3 & 4\\0 & 2\end{array}\right]
$$


<p>use</p>


```julia
    >>> Matrix([[1, -1], [3, 4], [0, 2]])
    ⎡1  -1⎤
    ⎢     ⎥
    ⎢3  4 ⎥
    ⎢     ⎥
    ⎣0  2 ⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>In <code>Julia</code>, the <code>Matrix</code> constructor is <em>not</em> exported, so must be qualified. Here we <em>avoid</em> thte vector of row vectors above:</p>
</li>
</ul>


```julia
sympy.Matrix([1 -1; 3 4; 0 2])
```




\[\left[ \begin{array}{rr}1&-1\\3&4\\0&2\end{array}\right]\]



<p><em>However</em>, through the magic of <code>PyCall</code>, such matrices are converted into <code>Julia</code> matrices, of type <code>Array&#123;Sym&#125;</code>, so the familiar matrix operations for <code>Julia</code> users are available.</p>

<p>In fact, the above could be done in the more <code>Julia</code>n manner through</p>


```julia
Sym[1 -1; 3 4; 0 2]
```




\[\left[ \begin{array}{rr}1&-1\\3&4\\0&2\end{array}\right]\]



<p>using an annotation to ensure the type. Alternatively, through promotion, just a single symbolic object will result in the same:</p>


```julia
[Sym(1) -1; 3 4; 0 2]
```




\[\left[ \begin{array}{rr}1&-1\\3&4\\0&2\end{array}\right]\]



<hr />

<p>To make it easy to make column vectors, a list of elements is considered to be a column vector.</p>


```julia
    >>> Matrix([1, 2, 3])
    ⎡1⎤
    ⎢ ⎥
    ⎢2⎥
    ⎢ ⎥
    ⎣3⎦
```

<h5>In <code>Julia</code>:</h5>

<p>For ths use, <code>sympy.Matrix</code> does work, but again its usage is discouraged:</p>


```julia
sympy.Matrix([1, 2, 3])
```




\[\left[ \begin{array}{r}1\\2\\3\end{array}\right]\]



<ul>
<li><p>Again, this is converted into a <code>Vector&#123;Sym&#125;</code> object or entered directly:</p>
</li>
</ul>


```julia
Sym[1,2,3]
```




\[ \left[ \begin{array}{r}1\\2\\3\end{array} \right] \]



<div class="admonition note"><p class="admonition-title">And again:</p><p>Using <code>sympy.Matrix</code> is strongly discouraged.</p>
</div>

<hr />

<p>Matrices are manipulated just like any other object in SymPy or Python.</p>


```julia
    >>> M = Matrix([[1, 2, 3], [3, 2, 1]])
    >>> N = Matrix([0, 1, 1])
    >>> M*N
    ⎡5⎤
    ⎢ ⎥
    ⎣3⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>In <code>Julia</code>, matrices are just matrices, and inherit all of the operations defined on them:</p>
</li>
</ul>


```julia
M = Sym[1 2 3; 3 2 1]
N = Sym[0, 1, 1]
M*N
```




    2-element Array{Any,1}:
     5
     3



<hr />

<p>One important thing to note about SymPy matrices is that, unlike every other object in SymPy, they are mutable.  This means that they can be modified in place, as we will see below.  The downside to this is that <code>Matrix</code> cannot be used in places that require immutability, such as inside other SymPy expressions or as keys to dictionaries.  If you need an immutable version of <code>Matrix</code>, use <code>ImmutableMatrix</code>.</p>

<h5>In <code>Julia</code>:</h5>

<p>A distinction is made between <code>ImmutableMatrix</code> and a mutable one. Mutable ones are mapped to <code>Julia</code> arrays, immutable ones are left as a symbolic object of type <code>SymMatrix</code>. The usual infix mathematical operations (but not dot broadcasting), 0-based indexing, and dot call syntax for methods are used with these objects.</p>

<h2>Basic Operations</h2>

<h3>Shape</h3>

<p>Here are some basic operations on <code>Matrix</code>.  To get the shape of a matrix use <code>shape</code></p>


```julia
    >>> M = Matrix([[1, 2, 3], [-2, 0, 4]])
    >>> M
    ⎡1   2  3⎤
    ⎢        ⎥
    ⎣-2  0  4⎦
    >>> M.shape
    (2, 3)
```

<h5>In <code>Julia</code>:</h5>


```julia
M = Sym[1 2 3; -2 0 4]
M
```




\[\left[ \begin{array}{rrr}1&2&3\\-2&0&4\end{array}\right]\]




```julia
M.shape
```




    (2, 3)



<p>Or, the <code>Julia</code>n counterpart:</p>


```julia
size(M)
```




    (2, 3)



<hr />

<h3>Accessing Rows and Columns</h3>

<p>To get an individual row or column of a matrix, use <code>row</code> or <code>col</code>.  For example, <code>M.row&#40;0&#41;</code> will get the first row. <code>M.col&#40;-1&#41;</code> will get the last column.</p>


```julia
    >>> M.row(0)
    [1  2  3]
    >>> M.col(-1)
    ⎡3⎤
    ⎢ ⎥
    ⎣4⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>these 0-based operations are supported:</p>
</li>
</ul>


```julia
M.row(0)
M.col(-1)
```




\[\left[ \begin{array}{r}3\\4\end{array}\right]\]



<p>The more familiar counterparts would be:</p>


```julia
M[1,:], M[:, end]
```




    (Sym[1, 2, 3], Sym[3, 4])



<hr />

<h3>Deleting and Inserting Rows and Columns</h3>

<p>To delete a row or column, use <code>row_del</code> or <code>col_del</code>.  These operations will modify the Matrix <strong>in place</strong>.</p>


```julia
    >>> M.col_del(0)
    >>> M
    ⎡2  3⎤
    ⎢    ⎥
    ⎣0  4⎦
    >>> M.row_del(1)
    >>> M
    [2  3]
```

<h5>In <code>Julia</code>:</h5>

<p>These methods do <strong>not</strong> work on <code>Array&#123;Sym&#125;</code> objects, use <code>Julia&#39;s</code> indexing notation to remove a row or column.</p>

<p>However, these methods <strong>do</strong> work on the <code>ImmutableMatrix</code> class:</p>


```julia
M = sympy.ImmutableMatrix([1 2 3; -2 0 4])  # avoid vector of row vector construction
M.col_del(0)
```




\begin{equation*}\left[\begin{matrix}2 & 3\\0 & 4\end{matrix}\right]\end{equation*}




```julia
M.row_del(1)
```




\begin{equation*}\left[\begin{matrix}1 & 2 & 3\end{matrix}\right]\end{equation*}



<p>!!! Alert</p>

<p>For older versions of sympy, the following did not work (using symbolic  values as matrix entries without reverting to  their PyObjects  had shape issues) but  this should work  now:</p>


```julia
@vars x
sympy.ImmutableMatrix([x 1;  1  x])
```




\begin{equation*}\left[\begin{matrix}x & 1\\1 & x\end{matrix}\right]\end{equation*}



<hr />

<div class="admonition note"><p class="admonition-title">TODO</p><p>This is a mess. See issue 6992.</p>
</div>

<p>To insert rows or columns, use <code>row_insert</code> or <code>col_insert</code>.  These operations <strong>do not</strong> operate in place.</p>


```julia
    >>> M
    [2  3]
    >>> M = M.row_insert(1, Matrix([[0, 4]]))
    >>> M
    ⎡2  3⎤
    ⎢    ⎥
    ⎣0  4⎦
    >>> M = M.col_insert(0, Matrix([1, -2]))
    >>> M
    ⎡1   2  3⎤
    ⎢        ⎥
    ⎣-2  0  4⎦
```

<h5>In <code>Julia</code>:</h5>


```julia
M
M = M.row_insert(1, Sym[0 4])
M
```




\begin{equation*}\left[\begin{matrix}1 & 2 & 3\\-2 & 0 & 4\end{matrix}\right]\end{equation*}




```julia
M = M.col_insert(0, Matrix([1, -2]))
M
```




\begin{equation*}\left[\begin{matrix}1 & 2 & 3\\-2 & 0 & 4\end{matrix}\right]\end{equation*}



<hr />

<p>Unless explicitly stated, the methods mentioned below do not operate in place. In general, a method that does not operate in place will return a new <code>Matrix</code> and a method that does operate in place will return <code>None</code>.</p>

<h5>In <code>Julia</code></h5>

<p>This would be the case for the immutable matrices.</p>

<hr />

<h2>Basic Methods</h2>

<p>As noted above, simple operations like addition and multiplication are done just by using <code>&#43;</code>, <code>*</code>, and <code>**</code>.  To find the inverse of a matrix, just raise it to the <code>-1</code> power.</p>


```julia
    >>> M = Matrix([[1, 3], [-2, 3]])
    >>> N = Matrix([[0, 3], [0, 7]])
    >>> M + N
    ⎡1   6 ⎤
    ⎢      ⎥
    ⎣-2  10⎦
    >>> M*N
    ⎡0  24⎤
    ⎢     ⎥
    ⎣0  15⎦
    >>> 3*M
    ⎡3   9⎤
    ⎢     ⎥
    ⎣-6  9⎦
    >>> M**2
    ⎡-5  12⎤
    ⎢      ⎥
    ⎣-8  3 ⎦
    >>> M**-1
    ⎡1/3  -1/3⎤
    ⎢         ⎥
    ⎣2/9  1/9 ⎦
    >>> N**-1
    Traceback (most recent call last):
    ...
    ValueError: Matrix det == 0; not invertible.
```

<h5>In <code>Julia</code>:</h5>


```julia
M = Sym[1 3; -2 3]
M1 = Sym[0 3; 0 7]
M + M1
```




\[\left[ \begin{array}{rr}1&6\\-2&10\end{array}\right]\]




```julia
M*M1
```




    2×2 Array{Any,2}:
     0  24
     0  15




```julia
3*M
```




\[\left[ \begin{array}{rr}3&9\\-6&9\end{array}\right]\]




```julia
M^2
```




    2×2 Array{Any,2}:
     -5  12
     -8   3




```julia
M^-1
```




\[\left[ \begin{array}{rr}\frac{1}{3}&- \frac{1}{3}\\\frac{2}{9}&\frac{1}{9}\end{array}\right]\]




```julia
M1^-1
```




    PyError ($(Expr(:escape, :(ccall(#= /Users/verzani/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'sympy.matrices.common.NonInvertibleMatrixError'>
    NonInvertibleMatrixError('Matrix det == 0; not invertible.')
      File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3672, in inv
        return self._eval_inverse(**kwargs)
      File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/dense.py", line 265, in _eval_inverse
        rv = M.inverse_GE(iszerofunc=iszerofunc)
      File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3587, in inverse_GE
        raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
    




<p>The above (except for the inverses) are using generic <code>Julia</code> definitions. For immutable matrices, we would have:</p>


```julia
M = sympy.ImmutableMatrix([1 3; -2 3])
M1 = sympy.ImmutableMatrix([0 3; 0 7])
M + M1
```




\begin{equation*}\left[\begin{matrix}1 & 6\\-2 & 10\end{matrix}\right]\end{equation*}




```julia
M*M1
```




\begin{equation*}\left[\begin{matrix}0 & 24\\0 & 15\end{matrix}\right]\end{equation*}




```julia
3*M
```




\begin{equation*}\left[\begin{matrix}3 & 9\\-6 & 9\end{matrix}\right]\end{equation*}




```julia
M^2
```




\begin{equation*}\left(\left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]\right)^{2}\end{equation*}




```julia
M^-1
```




\begin{equation*}\left[\begin{matrix}\frac{1}{3} & - \frac{1}{3}\\\frac{2}{9} & \frac{1}{9}\end{matrix}\right]\end{equation*}




```julia
M1^-1
```




    PyError ($(Expr(:escape, :(ccall(#= /Users/verzani/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'sympy.matrices.common.NonInvertibleMatrixError'>
    NonInvertibleMatrixError('Matrix det == 0; not invertible.')
      File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3672, in inv
        return self._eval_inverse(**kwargs)
      File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/dense.py", line 265, in _eval_inverse
        rv = M.inverse_GE(iszerofunc=iszerofunc)
      File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3587, in inverse_GE
        raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
    




<ul>
<li><p>There is no broadcasting defined for the <code>SymMatrix</code> type.</p>
</li>
</ul>

<hr />

<p>To take the transpose of a Matrix, use <code>T</code>.</p>


```julia
    >>> M = Matrix([[1, 2, 3], [4, 5, 6]])
    >>> M
    ⎡1  2  3⎤
    ⎢       ⎥
    ⎣4  5  6⎦
    >>> M.T
    ⎡1  4⎤
    ⎢    ⎥
    ⎢2  5⎥
    ⎢    ⎥
    ⎣3  6⎦
```

<h5>In <code>Julia</code>:</h5>


```julia
M = Sym[1 2 3; 4 5 6]
M
```




\[\left[ \begin{array}{rrr}1&2&3\\4&5&6\end{array}\right]\]




```julia
M.T
```




\[\left[ \begin{array}{rr}1&4\\2&5\\3&6\end{array}\right]\]



<hr />

<h2>Matrix Constructors</h2>

<p>Several constructors exist for creating common matrices.  To create an identity matrix, use <code>eye</code>.  The command <code>eye&#40;n&#41;</code> will create an <code>n x n</code> identity matrix:</p>


```julia
    >>> eye(3)
    ⎡1  0  0⎤
    ⎢       ⎥
    ⎢0  1  0⎥
    ⎢       ⎥
    ⎣0  0  1⎦
    >>> eye(4)
    ⎡1  0  0  0⎤
    ⎢          ⎥
    ⎢0  1  0  0⎥
    ⎢          ⎥
    ⎢0  0  1  0⎥
    ⎢          ⎥
    ⎣0  0  0  1⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p><code>eye</code> is <em>not</em> exported so must qualified:</p>
</li>
</ul>


```julia
sympy.eye(3)
sympy.eye(4)
```




\[\left[ \begin{array}{rrrr}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right]\]



<hr />

<p>To create a matrix of all zeros, use <code>zeros</code>.  <code>zeros&#40;n, m&#41;</code> creates an <code>n x m</code> matrix of <code>0</code>s.</p>


```julia
    >>> zeros(2, 3)
    ⎡0  0  0⎤
    ⎢       ⎥
    ⎣0  0  0⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>zeros is exported but the method expects a symbolic first argument.  Either qualify it:</p>
</li>
</ul>


```julia
sympy.zeros(2, 3)
```




\[\left[ \begin{array}{rrr}0&0&0\\0&0&0\end{array}\right]\]



<p><em>or</em> create a symbolic first value:</p>


```julia
zeros(Sym(2), 3)
```




\[\left[ \begin{array}{rrr}0&0&0\\0&0&0\end{array}\right]\]



<p><em>or</em> use the <code>Julia</code> constructor:</p>


```julia
zeros(Sym, 2, 3)
```




\[\left[ \begin{array}{rrr}0&0&0\\0&0&0\end{array}\right]\]



<hr />

<p>Similarly, <code>ones</code> creates a matrix of ones.</p>


```julia
    >>> ones(3, 2)
    ⎡1  1⎤
    ⎢    ⎥
    ⎢1  1⎥
    ⎢    ⎥
    ⎣1  1⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>Similarly with <code>ones</code>:</p>
</li>
</ul>


```julia
sympy.ones(3, 2)
```




\[\left[ \begin{array}{rr}1&1\\1&1\\1&1\end{array}\right]\]



<hr />

<p>To create diagonal matrices, use <code>diag</code>.  The arguments to <code>diag</code> can be either numbers or matrices.  A number is interpreted as a <code>1 x 1</code> matrix. The matrices are stacked diagonally.  The remaining elements are filled with <code>0</code>\ s.</p>


```julia
    >>> diag(1, 2, 3)
    ⎡1  0  0⎤
    ⎢       ⎥
    ⎢0  2  0⎥
    ⎢       ⎥
    ⎣0  0  3⎦
    >>> diag(-1, ones(2, 2), Matrix([5, 7, 5]))
    ⎡-1  0  0  0⎤
    ⎢           ⎥
    ⎢0   1  1  0⎥
    ⎢           ⎥
    ⎢0   1  1  0⎥
    ⎢           ⎥
    ⎢0   0  0  5⎥
    ⎢           ⎥
    ⎢0   0  0  7⎥
    ⎢           ⎥
    ⎣0   0  0  5⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>similarly with <code>diag</code>:</p>
</li>
</ul>


```julia
sympy.diag(1, 2, 3)
```




\[\left[ \begin{array}{rrr}1&0&0\\0&2&0\\0&0&3\end{array}\right]\]




```julia
sympy.diag(-1, sympy.ones(2, 2), sympy.Matrix([5, 7, 5]))
```




\[\left[ \begin{array}{rrrr}-1&0&0&0\\0&1&1&0\\0&1&1&0\\0&0&0&5\\0&0&0&7\\0&0&0&5\end{array}\right]\]



<ul>
<li><p>The first one, could also use <code>Julia</code>'s <code>diagm</code> function:</p>
</li>
</ul>


```julia
using LinearAlgebra
diagm(0 => Sym[1,2,3])
```




\[\left[ \begin{array}{rrr}1&0&0\\0&2&0\\0&0&3\end{array}\right]\]



<hr />

<h2>Advanced Methods</h2>

<h3>Determinant</h3>

<p>To compute the determinant of a matrix, use <code>det</code>.</p>


```julia
    >>> M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])
    >>> M
    ⎡1  0   1⎤
    ⎢        ⎥
    ⎢2  -1  3⎥
    ⎢        ⎥
    ⎣4  3   2⎦
    >>> M.det()
    -1
```

<h5>In <code>Julia</code>:</h5>


```julia
using LinearAlgebra
M = Sym[1 0 1; 2 -1 3; 4 3 2]
M
```




\[\left[ \begin{array}{rrr}1&0&1\\2&-1&3\\4&3&2\end{array}\right]\]




```julia
M.det()
```




\begin{equation*}-1\end{equation*}



<p>Let</p>


```julia
@vars x
A = Sym[x 1; 1 x]
```




\[\left[ \begin{array}{rr}x&1\\1&x\end{array}\right]\]



<p>Then we can compute the determinant using <code>Julia</code>'s generic implementation:</p>


```julia
det(A)
```




\begin{equation*}x^{2} - 1\end{equation*}



<p><em>or</em> using SymPy's:</p>


```julia
A.det()
```




\begin{equation*}x^{2} - 1\end{equation*}



<p>The answer is identical, though not necessarily being done in a similar manner.</p>

<h3>RREF</h3>

<p>To put a matrix into reduced row echelon form, use <code>rref</code>.  <code>rref</code> returns a tuple of two elements. The first is the reduced row echelon form, and the second is a tuple of indices of the pivot columns.</p>


```julia
    >>> M = Matrix([[1, 0, 1, 3], [2, 3, 4, 7], [-1, -3, -3, -4]])
    >>> M
    ⎡1   0   1   3 ⎤
    ⎢              ⎥
    ⎢2   3   4   7 ⎥
    ⎢              ⎥
    ⎣-1  -3  -3  -4⎦
    >>> M.rref()
    ⎛⎡1  0   1    3 ⎤        ⎞
    ⎜⎢              ⎥        ⎟
    ⎜⎢0  1  2/3  1/3⎥, (0, 1)⎟
    ⎜⎢              ⎥        ⎟
    ⎝⎣0  0   0    0 ⎦        ⎠
```

<h5>In <code>Julia</code>:</h5>


```julia
M = Sym[1 0 1 3; 2 3 4 7; -1 -3 -3 -4]
M
```




\[\left[ \begin{array}{rrrr}1&0&1&3\\2&3&4&7\\-1&-3&-3&-4\end{array}\right]\]




```julia
M.rref()
```




    (Sym[Sym(PyObject 1) Sym(PyObject 0) Sym(PyObject 1) Sym(PyObject 3); Sym(PyObject 0) Sym(PyObject 1) Sym(PyObject 2/3) Sym(PyObject 1/3); Sym(PyObject 0) Sym(PyObject 0) Sym(PyObject 0) Sym(PyObject 0)], (0, 1))



<div class="admonition note"><p class="admonition-title">Note</p><p>The first element of the tuple returned by <code>rref</code> is of type <code>Matrix</code>. The second is of type <code>tuple</code>.</p>
</div>

<h2>Nullspace</h2>

<p>To find the nullspace of a matrix, use <code>nullspace</code>. <code>nullspace</code> returns a <code>list</code> of column vectors that span the nullspace of the matrix.</p>


```julia
    >>> M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])
    >>> M
    ⎡1  2   3  0  0⎤
    ⎢              ⎥
    ⎣4  10  0  0  1⎦
    >>> M.nullspace()
    ⎡⎡-15⎤  ⎡0⎤  ⎡ 1  ⎤⎤
    ⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
    ⎢⎢ 6 ⎥  ⎢0⎥  ⎢-1/2⎥⎥
    ⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
    ⎢⎢ 1 ⎥, ⎢0⎥, ⎢ 0  ⎥⎥
    ⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
    ⎢⎢ 0 ⎥  ⎢1⎥  ⎢ 0  ⎥⎥
    ⎢⎢   ⎥  ⎢ ⎥  ⎢    ⎥⎥
    ⎣⎣ 0 ⎦  ⎣0⎦  ⎣ 1  ⎦⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>the list is mapped to an array of vectors, otherwise this is identical:</p>
</li>
</ul>


```julia
M = Sym[1 2 3 0 0; 4 10 0 0 1]
M
```




\[\left[ \begin{array}{rrrrr}1&2&3&0&0\\4&10&0&0&1\end{array}\right]\]




```julia
M.nullspace()
```




    3-element Array{Array{Sym,2},1}:
     [-15; 6; 1; 0; 0] 
     [0; 0; 0; 1; 0]   
     [1; -1/2; 0; 0; 1]



<h2>Columnspace</h2>

<p>To find the columnspace of a matrix, use <code>columnspace</code>. <code>columnspace</code> returns a <code>list</code> of column vectors that span the columnspace of the matrix.</p>


```julia
    >>> M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]])
    >>> M
    ⎡1  1  2⎤
    ⎢       ⎥
    ⎢2  1  3⎥
    ⎢       ⎥
    ⎣3  1  4⎦
    >>> M.columnspace()
    ⎡⎡1⎤  ⎡1⎤⎤
    ⎢⎢ ⎥  ⎢ ⎥⎥
    ⎢⎢2⎥, ⎢1⎥⎥
    ⎢⎢ ⎥  ⎢ ⎥⎥
    ⎣⎣3⎦  ⎣1⎦⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>as with <code>nullspace</code>, the return value is a vector of vectors:</p>
</li>
</ul>


```julia
M = Sym[1 1 2; 2 1 3; 3 1 4]
M
```




\[\left[ \begin{array}{rrr}1&1&2\\2&1&3\\3&1&4\end{array}\right]\]




```julia
M.columnspace()
```




    2-element Array{Array{Sym,2},1}:
     [1; 2; 3]
     [1; 1; 1]



<hr />

<h2>Eigenvalues, Eigenvectors, and Diagonalization</h2>

<p>To find the eigenvalues of a matrix, use <code>eigenvals</code>.  <code>eigenvals</code> returns a dictionary of <code>eigenvalue:algebraic multiplicity</code> pairs (similar to the output of :ref:<code>roots &lt;tutorial-roots&gt;</code>).</p>


```julia
    >>> M = Matrix([[3, -2,  4, -2], [5,  3, -3, -2], [5, -2,  2, -2], [5, -2, -3,  3]])
    >>> M
    ⎡3  -2  4   -2⎤
    ⎢             ⎥
    ⎢5  3   -3  -2⎥
    ⎢             ⎥
    ⎢5  -2  2   -2⎥
    ⎢             ⎥
    ⎣5  -2  -3  3 ⎦
    >>> M.eigenvals()
    {-2: 1, 3: 1, 5: 2}
```

<h5>In <code>Julia</code>:</h5>


```julia
M = Sym[3 -2  4 -2; 5  3 -3 -2; 5 -2  2 -2; 5 -2 -3  3]
M
```




\[\left[ \begin{array}{rrrr}3&-2&4&-2\\5&3&-3&-2\\5&-2&2&-2\\5&-2&-3&3\end{array}\right]\]




```julia
M.eigenvals()
```




    Dict{Any,Any} with 3 entries:
      3 => 1
      -2 => 1
      5 => 2



<hr />

<p>This means that <code>M</code> has eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.</p>

<p>To find the eigenvectors of a matrix, use <code>eigenvects</code>.  <code>eigenvects</code> returns a list of tuples of the form <code>&#40;eigenvalue:algebraic multiplicity, &#91;eigenvectors&#93;&#41;</code>.</p>


```julia
    >>> M.eigenvects()
    ⎡⎛       ⎡⎡0⎤⎤⎞  ⎛      ⎡⎡1⎤⎤⎞  ⎛      ⎡⎡1⎤  ⎡0 ⎤⎤⎞⎤
    ⎢⎜       ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥  ⎢  ⎥⎥⎟⎥
    ⎢⎜       ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥  ⎢-1⎥⎥⎟⎥
    ⎢⎜-2, 1, ⎢⎢ ⎥⎥⎟, ⎜3, 1, ⎢⎢ ⎥⎥⎟, ⎜5, 2, ⎢⎢ ⎥, ⎢  ⎥⎥⎟⎥
    ⎢⎜       ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥⎥⎟  ⎜      ⎢⎢1⎥  ⎢0 ⎥⎥⎟⎥
    ⎢⎜       ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥⎥⎟  ⎜      ⎢⎢ ⎥  ⎢  ⎥⎥⎟⎥
    ⎣⎝       ⎣⎣1⎦⎦⎠  ⎝      ⎣⎣1⎦⎦⎠  ⎝      ⎣⎣0⎦  ⎣1 ⎦⎦⎠⎦
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>the output is less than desirable, as there is no special <code>show</code> method</p>
</li>
<li><p>the <code>eigvals</code> and <code>eigvecs</code> methods present the output in the manner that <code>Julia</code>'s generic functions do:</p>
</li>
</ul>


```julia
M.eigenvects()
```




    3-element Array{Tuple{Sym,Int64,Array{Array{Sym,2},1}},1}:
     (-2, 1, [[0; 1; 1; 1]])              
     (3, 1, [[1; 1; 1; 1]])               
     (5, 2, [[1; 1; 1; 0], [0; -1; 0; 1]])



<p>compare with</p>


```julia
eigvecs(M)
```




\[\left[ \begin{array}{rrrr}0&1&1&0\\1&1&1&-1\\1&1&1&0\\1&1&0&1\end{array}\right]\]



<hr />

<p>This shows us that, for example, the eigenvalue 5 also has geometric multiplicity 2, because it has two eigenvectors.  Because the algebraic and geometric multiplicities are the same for all the eigenvalues, <code>M</code> is diagonalizable.</p>

<p>To diagonalize a matrix, use <code>diagonalize</code>. <code>diagonalize</code> returns a tuple <code>&#40;P, D&#41;</code>, where <code>D</code> is diagonal and <code>M &#61; PDP^&#123;-1&#125;</code>.</p>


```julia
    >>> P, D = M.diagonalize()
    >>> P
    ⎡0  1  1  0 ⎤
    ⎢           ⎥
    ⎢1  1  1  -1⎥
    ⎢           ⎥
    ⎢1  1  1  0 ⎥
    ⎢           ⎥
    ⎣1  1  0  1 ⎦
    >>> D
    ⎡-2  0  0  0⎤
    ⎢           ⎥
    ⎢0   3  0  0⎥
    ⎢           ⎥
    ⎢0   0  5  0⎥
    ⎢           ⎥
    ⎣0   0  0  5⎦
    >>> P*D*P**-1
    ⎡3  -2  4   -2⎤
    ⎢             ⎥
    ⎢5  3   -3  -2⎥
    ⎢             ⎥
    ⎢5  -2  2   -2⎥
    ⎢             ⎥
    ⎣5  -2  -3  3 ⎦
    >>> P*D*P**-1 == M
    True
```

<h5>In <code>Julia</code>:</h5>


```julia
P, D = M.diagonalize()
P
```




\[\left[ \begin{array}{rrrr}0&1&1&0\\1&1&1&-1\\1&1&1&0\\1&1&0&1\end{array}\right]\]




```julia
D
```




\[\left[ \begin{array}{rrrr}-2&0&0&0\\0&3&0&0\\0&0&5&0\\0&0&0&5\end{array}\right]\]




```julia
P*D*P^-1
```




    4×4 Array{Any,2}:
     3  -2   4  -2
     5   3  -3  -2
     5  -2   2  -2
     5  -2  -3   3




```julia
P*D*P^-1 == M
```




    true



<div class="admonition note"><p class="admonition-title">Quick Tip</p></div>

<p><code>lambda</code> is a reserved keyword in Python, so to create a Symbol called    $\lambda$, while using the same names for SymPy Symbols and Python    variables, use <code>lamda</code> (without the <code>b</code>).  It will still pretty print    as $\lambda$.</p>

<p>Note that since <code>eigenvects</code> also includes the eigenvalues, you should use it instead of <code>eigenvals</code> if you also want the eigenvectors. However, as computing the eigenvectors may often be costly, <code>eigenvals</code> should be preferred if you only wish to find the eigenvalues.</p>

<p>If all you want is the characteristic polynomial, use <code>charpoly</code>.  This is more efficient than <code>eigenvals</code>, because sometimes symbolic roots can be expensive to calculate.</p>


```julia
    >>> lamda = symbols('lamda')
    >>> p = M.charpoly(lamda)
    >>> factor(p)
           2
    (λ - 5) ⋅(λ - 3)⋅(λ + 2)
```

<h5>In <code>Julia</code>:</h5>

<ul>
<li><p>note missing <code>b</code> is not needed with <code>Julia</code>:</p>
</li>
</ul>


```julia
lambda = symbols("lambda")
p = M.charpoly(lambda)
factor(p)
```




\begin{equation*}\left(\lambda - 5\right)^{2} \left(\lambda - 3\right) \left(\lambda + 2\right)\end{equation*}



<hr />

<div class="admonition note"><p class="admonition-title">TODO</p><p>Add an example for <code>jordan_form</code>, once it is fully implemented.</p>
</div>

<h2>Possible Issues</h2>

<h2>Zero Testing</h2>

<p>If your matrix operations are failing or returning wrong answers, the common reasons would likely be from zero testing. If there is an expression not properly zero-tested, it can possibly bring issues in finding pivots for gaussian elimination, or deciding whether the matrix is inversible, or any high level functions which relies on the prior procedures.</p>

<p>Currently, the SymPy's default method of zero testing <code>_iszero</code> is only guaranteed to be accurate in some limited domain of numerics and symbols, and any complicated expressions beyond its decidability are treated as <code>None</code>, which behaves similarly to logical <code>False</code>.</p>

<p>The list of methods using zero testing procedures are as followings.</p>

<p><code>echelon_form</code> , <code>is_echelon</code> , <code>rank</code> , <code>rref</code> , <code>nullspace</code> , <code>eigenvects</code> , <code>inverse_ADJ</code> , <code>inverse_GE</code> , <code>inverse_LU</code> , <code>LUdecomposition</code> , <code>LUdecomposition_Simple</code> , <code>LUsolve</code></p>

<p>They have property <code>iszerofunc</code> opened up for user to specify zero testing method, which can accept any function with single input and boolean output, while being defaulted with <code>_iszero</code>.</p>

<p>Here is an example of solving an issue caused by undertested zero. [#zerotestexampleidea-fn]_ [#zerotestexamplediscovery-fn]_</p>


```julia
    >>> from sympy import *
    >>> q = Symbol("q", positive = True)
    >>> m = Matrix([
    ... [-2*cosh(q/3),      exp(-q),            1],
    ... [      exp(q), -2*cosh(q/3),            1],
    ... [           1,            1, -2*cosh(q/3)]])
    >>> m.nullspace()
    []
```

<h5>In <code>Julia</code>:</h5>


```julia
q = sympy.Symbol("q", positive = true)
m = Sym[
-2*cosh(q/3)      exp(-q)            1;
      exp(q) -2*cosh(q/3)            1;
           1            1 -2*cosh(q/3)]
m.nullspace()
```




    0-element Array{Any,1}



<hr />

<p>You can trace down which expression is being underevaluated, by injecting a custom zero test with warnings enabled.</p>


```julia
    >>> import warnings
    >>>
    >>> def my_iszero(x):
    ...     try:
    ...         result = x.is_zero
    ...     except AttributeError:
    ...         result = None
    ...
    ...     # Warnings if evaluated into None
    ...     if result == None:
    ...         warnings.warn("Zero testing of {} evaluated into {}".format(x, result))
    ...     return result
    ...
    >>> m.nullspace(iszerofunc=my_iszero) # doctest: +SKIP
    __main__:9: UserWarning: Zero testing of 4*cosh(q/3)**2 - 1 evaluated into None
    __main__:9: UserWarning: Zero testing of (-exp(q) - 2*cosh(q/3))*(-2*cosh(q/3) - exp(-q)) - (4*cosh(q/3)**2 - 1)**2 evaluated into None
    __main__:9: UserWarning: Zero testing of 2*exp(q)*cosh(q/3) - 16*cosh(q/3)**4 + 12*cosh(q/3)**2 + 2*exp(-q)*cosh(q/3) evaluated into None
    __main__:9: UserWarning: Zero testing of -(4*cosh(q/3)**2 - 1)*exp(-q) - 2*cosh(q/3) - exp(-q) evaluated into None
    []
```

<h5>In <code>Julia</code>:</h5>

<p>Is this available??</p>

<hr />

<p>In this case, <code>&#40;-exp&#40;q&#41; - 2*cosh&#40;q/3&#41;&#41;*&#40;-2*cosh&#40;q/3&#41; - exp&#40;-q&#41;&#41; - &#40;4*cosh&#40;q/3&#41;**2 - 1&#41;**2</code> should yield zero, but the zero testing had failed to catch. possibly meaning that a stronger zero test should be introduced. For this specific example, rewriting to exponentials and applying simplify would make zero test stronger for hyperbolics, while being harmless to other polynomials or transcendental functions.</p>


```julia
    >>> def my_iszero(x):
    ...     try:
    ...         result = x.rewrite(exp).simplify().is_zero
    ...     except AttributeError:
    ...         result = None
    ...
    ...     # Warnings if evaluated into None
    ...     if result == None:
    ...         warnings.warn("Zero testing of {} evaluated into {}".format(x, result))
    ...     return result
    ...
    >>> m.nullspace(iszerofunc=my_iszero) # doctest: +SKIP
    __main__:9: UserWarning: Zero testing of -2*cosh(q/3) - exp(-q) evaluated into None
    ⎡⎡  ⎛   q         ⎛q⎞⎞  -q         2⎛q⎞    ⎤⎤
    ⎢⎢- ⎜- ℯ  - 2⋅cosh⎜─⎟⎟⋅ℯ   + 4⋅cosh ⎜─⎟ - 1⎥⎥
    ⎢⎢  ⎝             ⎝3⎠⎠              ⎝3⎠    ⎥⎥
    ⎢⎢─────────────────────────────────────────⎥⎥
    ⎢⎢          ⎛      2⎛q⎞    ⎞     ⎛q⎞       ⎥⎥
    ⎢⎢        2⋅⎜4⋅cosh ⎜─⎟ - 1⎟⋅cosh⎜─⎟       ⎥⎥
    ⎢⎢          ⎝       ⎝3⎠    ⎠     ⎝3⎠       ⎥⎥
    ⎢⎢                                         ⎥⎥
    ⎢⎢           ⎛   q         ⎛q⎞⎞            ⎥⎥
    ⎢⎢          -⎜- ℯ  - 2⋅cosh⎜─⎟⎟            ⎥⎥
    ⎢⎢           ⎝             ⎝3⎠⎠            ⎥⎥
    ⎢⎢          ────────────────────           ⎥⎥
    ⎢⎢                   2⎛q⎞                  ⎥⎥
    ⎢⎢             4⋅cosh ⎜─⎟ - 1              ⎥⎥
    ⎢⎢                    ⎝3⎠                  ⎥⎥
    ⎢⎢                                         ⎥⎥
    ⎣⎣                    1                    ⎦⎦
```

<h5>In <code>Julia</code>:</h5>

<p>Is this available?</p>

<hr />

<p>You can clearly see <code>nullspace</code> returning proper result, after injecting an alternative zero test.</p>

<p>Note that this approach is only valid for some limited cases of matrices containing only numerics, hyperbolics, and exponentials. For other matrices, you should use different method opted for their domains.</p>

<p>Possible suggestions would be either taking advantage of rewriting and simplifying, with tradeoff of speed [#zerotestsimplifysolution-fn]_ , or using random numeric testing, with tradeoff of accuracy [#zerotestnumerictestsolution-fn]_ .</p>

<p>If you wonder why there is no generic algorithm for zero testing that can work with any symbolic entities, it's because of the constant problem stating that zero testing is undecidable [#constantproblemwikilink-fn]_ , and not only the SymPy, but also other computer algebra systems [#mathematicazero-fn]_ [#matlabzero-fn]_ would face the same fundamental issue.</p>

<p>However, discovery of any zero test failings can provide some good examples to improve SymPy, so if you have encountered one, you can report the issue to SymPy issue tracker [#sympyissues-fn]_ to get detailed help from the community.</p>

<div class="admonition note"><p class="admonition-title">Footnotes</p></div>

<ul>
<li><p>[#zerotestexampleidea-fn] Inspired by https://gitter.im/sympy/sympy?at=5b7c3e8ee5b40332abdb206c</p>
</li>
<li><p>[#zerotestexamplediscovery-fn] Discovered from https://github.com/sympy/sympy/issues/15141</p>
</li>
<li><p>[#zerotestsimplifysolution-fn] Suggested from https://github.com/sympy/sympy/issues/10120</p>
</li>
<li><p>[#zerotestnumerictestsolution-fn] Suggested from https://github.com/sympy/sympy/issues/10279</p>
</li>
<li><p>[#constantproblemwikilink-fn] https://en.wikipedia.org/wiki/Constant_problem</p>
</li>
<li><p>[#mathematicazero-fn] How mathematica tests zero https://reference.wolfram.com/language/ref/PossibleZeroQ.html</p>
</li>
<li><p>[#matlabzero-fn] How matlab tests zero https://www.mathworks.com/help/symbolic/mupad_ref/iszero.html</p>
</li>
<li><p>[#sympyissues-fn] https://github.com/sympy/sympy/issues</p>
</li>
</ul>

<hr />

<p><a href="./index.html">return to index</a></p>

1145ad9b09cf94ba7d4da4167100986fb4ff360e
ipynb
Jupyter Notebook
examples/matrices.ipynb
benneti/SymPy.jl
6a28bfd5286e6214ac3fec0f0ca2bdea4c0e1442
MIT
examples/matrices.ipynb
benneti/SymPy.jl
6a28bfd5286e6214ac3fec0f0ca2bdea4c0e1442
MIT
examples/matrices.ipynb
benneti/SymPy.jl
6a28bfd5286e6214ac3fec0f0ca2bdea4c0e1442
MIT
Qwen/Qwen-72B
1. YES
2. YES