413.sagews

I stayed up a bit and wrote an LU decomposition function for Sage.  You can see the actual code that is going into Sage soon at http://trac.sagemath.org/sage_trac/ticket/3048 or http://trac.sagemath.org/sage_trac/attachment/ticket/3048/lu_decomposition.patch.  In particular, you might notice that the new LU decomposition code is used in calculating the determinant of a matrix, solving a system of inequalities, and finding the inverse of a matrix.

%auto
%cython
from sage.matrix.matrix cimport Matrix
def _plu_decomposition_classical(Matrix self, pivoting=None):
    """
    Return the PLU decomposition P, L, and U such that self=P*L*U.
    If self is an mxn matrix, then P is an mxm permutation matrix,
    L is a mxm unit lower-triangular matrix, and U is an mxn
    upper-triangular matrix.

    The decomposition is done over the fraction field of self.base_ring().

    EXAMPLES::

    """
    #tm = verbose('generic in-place LU decomposition on %s x %s matrix'%(self._nrows, self._ncols))
    cdef Py_ssize_t start_row, c, r, nr, nc, i, pivot_row
    # Default to partial pivoting for inexact fields
    R = self.base_ring().fraction_field()

    x = self.fetch('PLU')
    if x is not None and x[3]==pivoting:
        return x[0].matrix(), x[1], x[2]

    cdef Matrix A, L

    nr = self._nrows
    nc = self._ncols
    R = self.base_ring().fraction_field()
    A = self.change_ring(R).copy()
    L = A.new_matrix(nr, nr)

    from sage.groups.all import SymmetricGroup 
    S = SymmetricGroup(nr)
    p = S(1)

    start_row = 0
    pivots = []
    for c from 0 <= c < nc:
        #print "column: ", c
        if PyErr_CheckSignals(): raise KeyboardInterrupt
        pivot = 0
        pivot_row = start_row
        if pivoting=="partial":
            # Pick largest possible pivot in column c
            for r from start_row <= r < nr:
                if pivot < abs(A.get_unsafe(r,c)):
                    pivot_row = r
                    pivot = abs(A.get_unsafe(pivot_row, c))
        else:
            # Pick the first nonzero pivot in column c
            for r from start_row <= r < nr:
                pivot = A.get_unsafe(r,c)
                if pivot != 0:
                    pivot_row = r
                    break

        if pivot != 0:
            pivots.append(c)
            if pivot_row != start_row:
                A.swap_rows(pivot_row, start_row)
                L.swap_rows(pivot_row, start_row)
                p *= S((pivot_row+1, start_row+1))
            a_inverse = ~A.get_unsafe(start_row,c)
            for i from start_row < i < nr:
                if A.get_unsafe(i,c):
                    b = A.get_unsafe(i, c)*a_inverse
                    A.add_multiple_of_row(i, start_row, -b, c)
                    L.set_unsafe(i,start_row,b)
            start_row += 1

    L += L.parent().identity_matrix()
    self.cache('pivots', pivots)
    self.cache('PLU', (p,L,A,pivoting))
    #verbose('done with LU decomposition', tm)            
    return p.matrix(), L, A

def lu_decomposition(Matrix self, *args, **kwds):
    """
    Wrapper function for PLU decomposition.
    """
    return _plu_decomposition_classical(self,*args,**kwds)


def solver(Matrix self, Matrix B, check_rank=True):
    r"""
    If self is a matrix `A` of full rank, then this function
    returns a matrix `X` such that `A X = B`.
    
    .. seealso::

       :meth:`solve_right` and :meth:`solve_left`
    
    INPUT:
    
    
    -  ``B`` - a matrix
    
    -  ``check_rank`` - bool (default: True) (currently ignored)
    
    
    OUTPUT: matrix
    
    EXAMPLES::
    
        sage: A = matrix(QQ,3,[1,2,4,5,3,1,1,2,-1])
        sage: B = matrix(QQ,3,2,[1,5,1,2,1,5])
        sage: A._solve_right_nonsingular_square(B)
        [ -1/7 -11/7]
        [  4/7  23/7]
        [    0     0]
        sage: A._solve_right_nonsingular_square(B, check_rank=False)
        [ -1/7 -11/7]
        [  4/7  23/7]
        [    0     0]
        sage: X = A._solve_right_nonsingular_square(B, check_rank=False)
        sage: A*X == B
        True
    """
    cdef Matrix L,U
    cdef int i,j
    # Should use LU Decomposition!
    lu_decomposition(self)
    P,L,U,strategy=self.fetch('PLU')
    # Swap matrix rows of B using P
    B=B[(~P)(range(B._nrows))].copy()
    # Forward substitution
    for i from 1<=i<B._nrows:
        B[i] = B[i]-sum([L[i,j]*B[j] for j in range(i)])
    # Backwards substitution
    for i from B._nrows>i>=0:
        B[i] = (B[i]-sum([U[i,j]*B[j] for j in range(U._ncols-1,i,-1)]))/U[i,i]
    return B
    
    
def det(Matrix self):
    cdef Matrix u
    lu_decomposition(self) 
    p,l,u,strategy=self.fetch('PLU') 
    d = p.sign() 
    for i from 0<= i < self._nrows: 
        d *= u.get_unsafe(i,i) 
    self.cache('det', d) 
    return d

Let's try the example from the book.

A=matrix(4,5,[2,4,-1,5,-2, -4,-5,3,-8,1, 2, -5, -4, 1, 8, -6, 0, 7, -3, 1])
show(A)

\left($$\begin{array}{rrrrr} 2 & 4 & -1 & 5 & -2 \\ -4 & -5 & 3 & -8 & 1 \\ 2 & -5 & -4 & 1 & 8 \\ -6 & 0 & 7 & -3 & 1 \end{array}$$\right)

If the following command doesn't work, evaluate the first cell.  You can also look at the source code in that cell for the LU decomposition.  It's slightly more general than the algorithm we talked about in class, but the underlying principles are exactly the same.  The generalization came in when doing the pivoting and row swaps.

P,L,U=lu_decomposition(A)
show(P)
show(L)
show(U)

\left($$\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array}$$\right)

\left($$\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 1 & -3 & 1 & 0 \\ -3 & 4 & 2 & 1 \end{array}$$\right)

\left($$\begin{array}{rrrrr} 2 & 4 & -1 & 5 & -2 \\ 0 & 3 & 1 & 2 & -3 \\ 0 & 0 & 0 & 2 & 1 \\ 0 & 0 & 0 & 0 & 5 \end{array}$$\right)

Let's check that it is right.

P*L*U==A

True

Remember that Matlab gave us a different answer?  That was because Matlab used partial pivoting (it chose the largest pivot and did a row swap).  Partial pivoting helps minimize numerical errors.  In Sage, because there is exact arithmetic (and we are dealing with exact numbers), we don't need to use partial pivoting.  However, we can still ask for it anyway.

P,L,U=lu_decomposition(A, pivoting="partial")
show(P)
show(L)
show(U)

\left($$\begin{array}{rrrr} 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \end{array}$$\right)

\left($$\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ \frac{2}{3} & 1 & 0 & 0 \\ -\frac{1}{3} & 1 & 1 & 0 \\ -\frac{1}{3} & -\frac{4}{5} & -\frac{2}{15} & 1 \end{array}$$\right)

\left($$\begin{array}{rrrrr} -6 & 0 & 7 & -3 & 1 \\ 0 & -5 & -\frac{5}{3} & -6 & \frac{1}{3} \\ 0 & 0 & 0 & 6 & 8 \\ 0 & 0 & 0 & 0 & -\frac{1}{3} \end{array}$$\right)

Notice that from $P$, we see that the first and last rows are swapped in order to get better pivots numerically.  Let's still verify the results.

P*L*U==A

True

You might notice in other examples that the Matlab $P$ is still a bit different than the Sage $P$.  That's because the Matlab $L$, $U$, and $P$ satisfy $LU=PA$, while the Sage results satisfy $PLU=A$.  As we'll find out later in the course, basically the Matlab $P$ is the transpose of the Sage $P$.  There are different conventions for where $P$ should go.  Some other software does a $LUP=A$ decomposition.  The Lapack software does a $PLU=A$ decomposition, like Sage.

$CC$ is the set of complex numbers, so here, we construct a random $4\times 5$ complex matrix.

a=random_matrix(CC,4,5)
show(a)

\left($$\begin{array}{rrrrr} 0.229692748016513 - 0.845661923562603i & 0.287536701715982 + 0.236946603380285i & -0.185534572565047 - 0.289195161295678i & 0.381784037827768 - 0.869211099061085i & -0.0586247227789061 - 0.293832522677727i \\ 0.572588741325297 - 0.188787845429612i & -0.179735519225174 + 0.226375704374467i & 0.807922936521494 + 0.513268974988697i & -0.652204136019794 + 0.672635212571936i & -0.0441922790787861 - 0.0549952757621675i \\ -0.630452932491459 + 0.132255191298339i & -0.154287888785182 + 0.924044425657864i & -0.533770069467662 + 0.866141745963207i & -0.341068398547977 - 0.107524013242364i & 0.722668222090292 - 0.299564064600013i \\ -0.156088117010016 + 0.654210787293331i & 0.194396079337710 + 0.910673000800208i & 0.652747709221585 + 0.354757282827718i & -0.415673580873110 + 0.342933310778510i & -0.0176146945731379 - 0.323775635584187i \end{array}$$\right)

If we have a matrix that is not exact (like this floating point complex matrix), the function above knows automatically to do partial pivoting to minimize numerical error.

P,L,U=lu_decomposition(a)
show(P)
show(L)
show(U)

\left($$\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{array}$$\right)

\left($$\begin{array}{rrrr} 1.00000000000000 & 0 & 0 & 0 \\ -0.334226714856104 - 0.654733841584202i & 1.00000000000000 & 0 & 0 \\ 0.379175796722956 + 0.574100355079673i & -0.000975088882645078 + 0.128358851883204i & 1.00000000000000 & 0 \\ -0.767145822259989 + 0.0237916766424910i & 0.821587531362721 - 0.500109006499683i & -0.0481482063416919 - 0.705230791651185i & 1.00000000000000 \end{array}$$\right)

\left($$\begin{array}{rrrrr} 0.229692748016513 - 0.845661923562603i & 0.287536701715982 + 0.236946603380285i & -0.185534572565047 - 0.289195161295678i & 0.381784037827768 - 0.869211099061085i & -0.0586247227789061 - 0.293832522677727i \\ 0 & -0.213322401451592 + 1.19149831981292i & -0.406434821225628 + 0.648009233808876i & 0.355635948235587 - 0.148070653655511i & 0.895456369941713 - 0.436154333329313i \\ -5.55111512312578 \times 10^{-17}i & -2.77555756156289 \times 10^{-17} - 3.46944695195361 \times 10^{-18}i & 0.795027522292052 + 0.782241618288173i & -1.31464120579083 + 0.737243267846046i & -0.245763678633482 - 0.0252896615076734i \\ 5.55111512312578 \times 10^{-17} & 5.55111512312578 \times 10^{-17} + 2.22044604925031 \times 10^{-16}i & 0 & -0.944828403398494 - 0.925080408664571i & -0.581148284506794 + 0.0838337447132584i \end{array}$$\right)

Check the results.  We can see that the difference between the two matrices is really small.

P*L*U-a

[ 0 0 0 0 0] [ 0 -2.77555756156289e-17 0 0 0] [ 0 0 0 0 0] [ -2.77555756156289e-17 0 0 -1.11022302462516e-16*I 0]

def make_symbolic_matrix(size):
    return matrix(size,[var('x%s%s'%(i,j)) for i in [0..size-1] for j in [0..size-1]])

A=make_symbolic_matrix(5)
A

[x00 x01 x02 x03 x04] [x10 x11 x12 x13 x14] [x20 x21 x22 x23 x24] [x30 x31 x32 x33 x34] [x40 x41 x42 x43 x44]

%time
c=A.det()

CPU time: 0.00 s, Wall time: 0.00 s

%time
b=det(A)

CPU time: 0.05 s, Wall time: 0.05 s

%time
b.expand()