SymPy: Extract the lower triangular part of a matrix

Question

I am trying to extract the lower triangular part of a SymPy matrix. Since I could not find a tril method in SymPy, I defined:

def tril (M):
    m = M.copy()
    for row_index in range (m.rows):
        for col_index in range (row_index + 1, m.cols):
            m[row_index, col_index] = 0
    return (m)

It seems to work:

Is there a more elegant way to extract the lower triangular part of a SymPy matrix?

Is .copy() the recommended way to ensure the integrity of the original matrix?

Answer 1

In SymPy, M.lower_triangular(k) will give the lower triangular elements below the kth diagonal. The default is k=0.

Answer 2

In [99]: M
Out[99]: 
?a  b  c?
?       ?
?d  e  f?
?       ?
?g  h  i?

The other answer suggest using the np.tril function:

In [100]: np.tril(M)
Out[100]: 
array([[a, 0, 0],
       [d, e, 0],
       [g, h, i]], dtype=object)

That converts M into a numpy array - object dtype because of the symbols. And the result is also a numpy array.

Your function returns a sympy.Matrix.

In [101]: def tril (M):
     ...:     m = M.copy()
     ...:     for row_index in range (m.rows):
     ...:         for col_index in range (row_index + 1, m.cols):
     ...:             m[row_index, col_index] = 0
     ...:     return (m)
     ...: 

In [102]: tril(M)
Out[102]: 
?a  0  0?
?       ?
?d  e  0?
?       ?
?g  h  i?

As a general rule mixing sympy and numpy leads to confusion, if not errors. numpy is best for numeric work. It can handle non-numeric objects like symbols, but the math is hit-or-miss.

The np.tri... functions are built on the np.tri function:

In [114]: np.tri(3).astype(int)
Out[114]: 
array([[1, 0, 0],
       [1, 1, 0],
       [1, 1, 1]])

We can make a symbolic Matrix from this:

In [115]: m1 = Matrix(np.tri(3).astype(int))

In [116]: m1
Out[116]: 
?1  0  0?
?       ?
?1  1  0?
?       ?
?1  1  1?

and do element-wise multiplication:

In [117]: M.multiply_elementwise(m1)
Out[117]: 
?a  0  0?
?       ?
?d  e  0?
?       ?
?g  h  i?

np.tri works by comparing a column array with a row:

In [123]: np.arange(3)[:,None]>=np.arange(3)
Out[123]: 
array([[ True, False, False],
       [ True,  True, False],
       [ True,  True,  True]])

In [124]: _.astype(int)
Out[124]: 
array([[1, 0, 0],
       [1, 1, 0],
       [1, 1, 1]])

Another answer suggests lower_triangular. It's interesting to look at its code:

    def entry(i, j):
        return self[i, j] if i + k >= j else self.zero

    return self._new(self.rows, self.cols, entry)

It applies an i>=j test to each element. _new must be iterating on the rows and columns.

Answer 3

You can simply use numpy function:

import numpy as np    
np.tril(M)

*of course, as noted below, you should convert back to sympy.Matrix(np.tril(M)). But it depends on what you're going to do next.