iconEuler Home

Linear Algebra

This is an introduction to the way EMT handles matrices for linear algebra. You should be familiar with the Euler matrix language by now.

01 - The Matrix Language

Let us start with a simple matrix.

>A:=[1,2,3;4,5,6;7,8,9]
            1             2             3 
            4             5             6 
            7             8             9 

det() computes the determinant, which surprisingly is 0 for this matrix.

>det(A)
0

Consequently, the kernel of A will contain non-zeros elements.

>kernel(A)
            1 
           -2 
            1 

Let us check, using the dot product.

>A.kernel(A)
            0 
            0 
            0 

By the way, the matrix language provides a very simple way to produce this matrix. We generate the vector 1:9, and put its entries into the rows of a 3x3 matrix.

>redim(1:9,3,3)
            1             2             3 
            4             5             6 
            7             8             9 

Now we try to solve a linear system.

>b:=[3;4;3]
            3 
            4 
            3 

This fails, since the matrix is singular.

>x:=A\b
Determinant zero!
Error in:
x:=A\b ...
      ^

Singularity is checked with the internal epsilon. In this case, the numerical solver recognizes that the matrix is singular.

EMT has another solver, which applies the pseudo-inverse of A to b. The solution is the best possible in distance with least norm. But it is not a proper solution.

>x:=svdsolve(A,b); fracprint(x)
     -5/3 
        0 
      5/3 

The error is still substantial, but it is the smallest possible error, nevertheless.

>norm(A.x-b)
0.816496580928

This is as close as we can get to b.

>A.x
      3.33333 
      3.33333 
      3.33333 

The fit function finds another solution with the same error. It is a bit faster.

>x:=fit(A,b); fracprint(x), norm(A.x-b)
    -10/3 
     10/3 
        0 
0.816496580928

We change one element of A.

>A[3,3]:=0
            1             2             3 
            4             5             6 
            7             8             0 

Now, A is regular.

>det(A)
27

And we can properly solve the linear system.

>x:=A\b; fraction x
    -19/9 
     20/9 
      2/9 

Indeed, this is the solution.

>norm(A.x-b)
0

Let us try a larger problem. We generate a 100x100 random matrix.

>A:=normal(100,100);

The determinant is huge.

>det(A)
-1.80259873511e+078

We multiply with a random vector.

>x:=random(rows(A),1); b=A.x;

The norm of the error is usually quite good. We need to switch to longestformat to avoid the error being rounded to 0 for the output.

>longest norm(A\b-x)
 8.661666511000772e-014 

Orthogonal Projection

In the following commands we demonstrate how to do the Gram-Schmid method by hand for the vectors v and w.

>v:=[1;1;1]; w:=[1;2;-1]; v1:=w-(w'.v)/(v'.v)*v; fracprint(v1);
      1/3 
      4/3 
     -5/3 

We could continue like this. But we take the cross product to get the third vector.

>v2:=crossproduct(v,w),
           -3 
            2 
            1 

Now, the normalized vectors v, v1, v2 should form an orthogonal matrix.

>M:=[v/norm(v),v1/norm(v1),v2/norm(v2)]
      0.57735      0.154303     -0.801784 
      0.57735      0.617213      0.534522 
      0.57735     -0.771517      0.267261 

Let us check this.

>norm(M.M'-id(3))
0

We can directly expand v to an orthonormal basis. For this, we note the space kernel(A') is orthogonal to all columns of a matrix A. Then we apply the utility function orthogonal() to all vectors.

>orthogonal(v|kernel(v'))
      0.57735     -0.707107     -0.408248 
      0.57735      0.707107     -0.408248 
      0.57735             0      0.816497 

Let us project the vector b to the plane generated by v and w. The coefficients of v and w minimize l1*v1+l2*v2-b.

The function "fit" does just this projection.

>A:=v|w; b:=[1;0;0]; l:=fit(A,b); fracprint(l);
      2/7 
     1/14 

We can compute the projection A.l.

>z:=A.l; fracprint(z);
     5/14 
      3/7 
     3/14 

And we compute the distance of z to this plane.

>norm(z-b)
0.801783725737

We could also use our matrix M to get the same result. The first two columns span the plane and are orthonormal.

>l:=M'.b; l=l[[1,2]]; M[:,[1,2]].l
     0.357143 
     0.428571 
     0.214286 

Another idea is to use the projection to v2. The matrix of this projection is v2.v2', after v2 is normed to 1.

>P:=id(3)-v2.v2'/norm(v2)^2; P.b
     0.357143 
     0.428571 
     0.214286 

Symbolic Algebra

Maxima knows a lot of functions for linear Algebra too. We can define the Matrix as data for symbolic expressions.

>A &:= [1,2,3;4,5,6;7,8,9]
            1             2             3 
            4             5             6 
            7             8             9 

The rank of the matrix is 2, because it is singular.

>&rank(A)
                                  2

The result of the Gauss algorithm applied to A can be computed with echelon.

The echelon form can be used to determine a base of the lines of A taking the non-zero lines of A. Moreover, the independent columns of the echelon form, are the columns of A, which form a basis of the columns of A.

>&echelon(A)
                             [ 1  2  3 ]
                             [         ]
                             [ 0  1  2 ]
                             [         ]
                             [ 0  0  0 ]

Here is one way to make a linear system from the matrix A, and to solve it. The system is A.[x,y,z]=[1,1,1].

>&transpose(A.[x,y,z]-[1,1,1]), sol &= solve(%[1],[x,y,z])
       [ 3 z + 2 y + x - 1  6 z + 5 y + 4 x - 1  9 z + 8 y + 7 x - 1 ]


               [[x = %r1 - 1, y = 1 - 2 %r1, z = %r1]]

The list of symbols is contained in a temporary list named %rnum_list, which we can access with temp(1). We can assign a value for each symbol as follows.

>&sol with temp(1)=0
                      [[x = - 1, y = 1, z = 0]]

For an invertible matrix, there is also a function lusolve(), which directly solves Ax=b.

>&lusolve([1,2,3;4,5,6;7,9,10],[1;1;1])
                               [  0  ]
                               [     ]
                               [ - 1 ]
                               [     ]
                               [  1  ]

With Maxima, we can use symbolic matrices containing variables or other expressions.

>A &= [1,2,3;4,5,6;7,8,a]
                             [ 1  2  3 ]
                             [         ]
                             [ 4  5  6 ]
                             [         ]
                             [ 7  8  a ]

We can easily see, that a=9 makes a singular matrix.

>&determinant(A)|factor
                             - 3 (a - 9)

We can now invert this matrix symbolically. This works unless the determinant is 0.

>Ainv &= ratsimp(inv(A)*det(A))
                     [ 5 a - 48  24 - 2 a  - 3 ]
                     [                         ]
                     [ 42 - 4 a   a - 21    6  ]
                     [                         ]
                     [   - 3        6      - 3 ]

We can print this with latex using the comment line "maxima: Ainv"

04 - Linear Algebra

Eigenvalues

Let us do some eigenvalue problems.

First we generate a diagonal matrix of size 3x3, with the specified values on the 0-th diagonal.

>D:=diagmatrix([1,2,3])
            1             0             0 
            0             2             0 
            0             0             3 

We generate a random matrix similar to the matrix A.

>M:=random(3,3); A:=M.D.inv(M)
      5.61971       3.54267      -5.38696 
      5.63576       6.88812      -8.29591 
      5.76664       5.13393      -6.50783 

Indeed, the eigenvalues of A are just 1,2,3. However, we get complex eigenvalues, which need to make real.

>real(eigenvalues(A))
[3,  1,  2]

We can also use the function eigen() to compute both the eigenvalues and the eigenspaces.

>{l,v}:=eigen(A); real(v),
    -0.863456      0.292455             1 
    -0.882092             1      0.108208 
           -1      0.908439      0.743101 

Then inv(v).A.v should be the diagonal matrix of eigenvalues.

>norm(inv(v).A.v-diagmatrix(l))
0

Euler can also do singular value decompositions.

>{p,s,q}:=svd(A); s
[17.8852,  1.79654,  0.186732]

These are the eigenvalues of A'.A, but Euler uses a special algorithm for this.

>sqrt(abs(eigenvalues(A'.A)))
[17.8852,  1.79654,  0.186732]

We get A as a composition of two orthogonal matrices, and a diagonal matrix.

>norm(p.diagmatrix(s).q'-A)
0

Matrix Manipulation

Let us demonstrate, how to solve a linear system by hand. Assume, we have the following system.

>fracformat(10); A:=redim(1:12,3,4)
        1         2         3         4 
        5         6         7         8 
        9        10        11        12 
>pivotize(A,1,1)
        1         2         3         4 
        0        -4        -8       -12 
        0        -8       -16       -24 
>pivotize(A,3,2)
        1         0        -1        -2 
        0         0         0         0 
        0         1         2         3 
>swapRows(A,2,3)
        1         0        -1        -2 
        0         1         2         3 
        0         0         0         0 

There is the function echelon, which does all this automatically.

>fracformat(10); A:=redim(1:12,3,4); echelon(A)
        1         0        -1        -2 
        0         1         2         3 

Eigenvalues with Maxima

We can also compute eigenvalues with Maxima. The function returns the eigenvalues and their multiplicities.

>longformat; A&:=[1,2;3,4]
                  1                   2 
                  3                   4 
>&eigenvalues(A)
                  5 - sqrt(33)  sqrt(33) + 5
                [[------------, ------------], [1, 1]]
                       2             2

Here are the eigenvectors.

>&eigenvectors(A)
          5 - sqrt(33)  sqrt(33) + 5
       [[[------------, ------------], [1, 1]], 
               2             2
                                 3 - sqrt(33)         sqrt(33) + 3
                           [[[1, ------------]], [[1, ------------]]]]
                                      4                    4

Large and Sparse Systems

Euler has a special data types, and algorithms for large and sparse system. For more information refer to the following introduction.

19 - Large Systems

The following does not use these algorithms, and the matrix is full. Yet, Euler can produce a good solution within seconds.

>A=normal(500,500); b=sum(A); x=A\b;
>totalmax(A.x-b)
0

Euler Home