# Linear algebra in Darwin

A vector in Darwin is delimited by square brackets: [....]

`v1 := [2,3,5,7,11,13];`
```v1 := [2, 3, 5, 7, 11, 13]
```

Arithmetic with vectors works as expected, provided that the lengths of the vectors are consistent. More precisely, addition and subtraction of vectors is done component-wise, provided that the vectors are of the same length and contain only numbers.

`v2 := 3*v1 - [5,5,5,5,5,5];`
```v2 := [1, 4, 10, 16, 28, 34]
```
`v1 + [1,1,1];`
```Error, different length arrays in sum
```

Multiplication or division of vectors by a scalar (a number) is done component-wise.

`v2/10;  100*v1;`
```[0.10000000, 0.4000, 1, 1.6000, 2.8000, 3.4000]
[200, 300, 500, 700, 1100, 1300]
```

Multiplication of vectors by vectors is understood to be an inner product and the vectors have to be of the same length. Vectors are both (simultaneously) column or row. The decision is made in the context of the operation.

```v1*v2;
(v1*v2) / sqrt( v1^2 * v2^2 );```
```926
0.9916
```

Notice that a vector squared (to the power 2) is equivalent to the inner product of the vector with itself, and hence equal to the norm squared.

A matrix is a list of vectors (of the same length).

```A := [ [2,3,4],
[-1,0,2],
[3,0,-1] ]:
print( A );```
```  2  3  4
-1  0  2
3  0 -1
```

The list components are the rows of the matrix, that is selecting a matrix with a single index will yield the corresponding row.

`A;`
```[-1, 0, 2]
```

Sometimes we want to generate a vectors or matrices with large dimensions, or dimensions which are given by some other variable(s). In this case, it is not suitable to build them with lists and the CreateArray function is used. CreateArray builds a new array and initializes all the elements to 0 (or to some other specified value)

`B := CreateArray(1..7);`
```B := [0, 0, 0, 0, 0, 0, 0]
```
`B := CreateArray( 1..3, 1..4, 77 ):  print(B);`
``` 77 77 77 77
77 77 77 77
77 77 77 77
```

Quite often it is very useful to build arrays with random values in them.

`C := Rand( array(integer,3,4) ):  print(C);`
``` -13  -6   6   1
0   0 -19   2
-3  -3  -2   1
```
`Rand( array(numeric,8) );`
```[0.5612, 0.5340, 0.1464, 0.08147475, 0.1956, 0.5040, 0.4691, 0.07674076]
```

Arithmetic with matrices works as with vectors. The normal display of matrices is as a list of list, which is not easy to read. By printing the matrix with `print` matrices are easier to read.

`print( 100*A );`
```  200  300  400
-100    0  200
300    0 -100
```
`A^2;  A^3;`
```[[13, 6, 10], [4, -3, -6], [3, 9, 13]]
[[50, 39, 54], [-7, 12, 16], [36, 9, 17]]
```
`print( A^3-A^2 );`
```  37  33  44
-11  15  22
33   0   4
```

A matrix can be multiplied by a vector on the left or on the right. When the vector is on the right, it is assumed that it is a column vector. When it is on the left it is assumed to be a row vector. The corresponding dimensions have to match.

```A*[1,1,1];
[1,1,1]*A;```
```[9, 1, 2]
[4, 3, 5]
```
`A*[1,2,3,4];`
```Error, unequal length vectors for inner product
```

The Fibonacci numbers are generated by the powers of the following matrix:

`F := [ [1,1], [1,0] ]:  print(F);`
``` 1 1
1 0
```
`for i to 5 do lprint( i, F^i ) od;`
```1 [[1, 1], [1, 0]]
2 [[2, 1], [1, 1]]
3 [[3, 2], [2, 1]]
4 [[5, 3], [3, 2]]
5 [[8, 5], [5, 3]]
```

notice that the successive powers of the matrix F contain the Fibonacci sequence.

Matrix inverses and transpositions are indicated by exponentiation. Matrix inverses can also be computed as division. Inverse computation can only be done on square, non-singular matrices. Transposition is indicated by the exponent t or T.

`Ainv := A^(-1):  print( Ainv );`
```  0.00000000  0.20000000  0.40000000
0.33333333 -0.93333333 -0.53333333
0.00000000  0.60000000  0.20000000
```
`print( A * Ainv );`
```  1.00000000 -0.00000000  0.00000000
0.00000000  1.00000000 -0.00000000
0.00000000 -0.00000000  1.00000000
```
`print( A^t );`
```  2 -1  3
3  0  0
4  2 -1
```

Notice that if a matrix is singular, computing its inverse will give an error:

```S := [ [1,2], [3,6] ]:
print(S);
1/S;```
``` 1 2
3 6
Error, (in list_power) matrix_inverse: matrix is singular
```

Gaussian elimination is a method to solve a linear system of equations. Given a matrix `A` and a vector `b`, gaussian elimination finds a vector `x` such that `Ax=b`. This can be done with a function call

```x := GaussElim( A, [1,2,3] );
A * x;```
```x := [1.6000, -3.1333, 1.8000]
[1.0000, 2.0000, 3]
```

or by computing the inverse and multiplying. Notice that the order of multiplication matters and the inverse should be multiplied on the left.

`1/A * [1,2,3];`
```[1.6000, -3.1333, 1.8000]
```

Symmetric matrices are very common and Darwin can compute eigenvalues/eigenvectors for symmetric matrices. First we will build a symmetric matrix with random integers.

```S := CreateArray( 1..4, 1..4 ):
for i to 4 do
for j from i to 4 do
S[i,j] := S[j,i] := Rand(integer) od od;
print(S);```
``` -24  -8  -1  -7
-8 -12   6  -5
-1   6  16  -2
-7  -5  -2  -7
```

The Eigenvalues function computes the eigenvalues and optionally the eigenvectors. It is easier to read its help file:

`?Eigenvalues`
```
Function Eigenvalues - Eigenvalue/vector decomposition of a symmetric matrix
Option: builtin

Calling Sequence:  Eigenvalues(A,eigenvects)

Parameters:
Name         Type     Description
---------------------------------------------
A            matrix   a symmetic matrix
eigenvects   name     an optional matrix name

Returns:
list(numeric)
Synopsis:  Compute an eigenvalue/eigenvector decomposition of A.  A must be a
symmetric matrix.  The function returns the vector containing the
eigenvalues in increasing order.  The optional second argument, if present
must be a name that will be assigned with the matrix of the eigenvectors.
The eigenvectors have norm 1 and are stored columnwise and the ith column
corresponds to the ith eigenvalue.

Examples:
> A := [[3,1,2],[1,2,-1],[2,-1,5]];
A := [[3, 1, 2], [1, 2, -1], [2, -1, 5]]
> alpha := Eigenvalues(A,V);
alpha := [0.4921, 3.2444, 6.2635]
> Vt := V^t;
Vt := [[0.6041, -0.6782, -0.4185], [0.6191, 0.7301, -0.2894], [0.5018, -0.08423029, 0.8609]]
> A*Vt = alpha*Vt;
[0.2973, -0.3337, -0.2059] = [0.2973, -0.3337, -0.2059]
> Vt*Vt;
1.0000
?Cholesky    ?GaussElim     ?Identity    ?matrix_inverse    ?transpose
?convolve    ?GivensElim    ?matrix      ?SvdAnalysis

------------
```

We can now compute the eigenvalues and eigenvectors and verify the main theorem of the eigenvalue/eigenvector decomposition.

```lam := Eigenvalues( S, U );
LAM := CreateArray( 1..4, 1..4 ):
for i to 4 do LAM[i,i] := lam[i] od:
print(LAM);```
```lam := [-31.0442, -9.4038, -4.2395, 17.6876]
-31.0442466   0.0000000   0.0000000   0.0000000
0.0000000  -9.4038345   0.0000000   0.0000000
0.0000000   0.0000000  -4.2394718   0.0000000
0.0000000   0.0000000   0.0000000  17.6875529
```
`print( U * LAM * U^t );`
``` -24.0000000  -8.0000000  -1.0000000  -7.0000000
-8.0000000 -12.0000000   6.0000000  -5.0000000
-1.0000000   6.0000000  16.0000000  -2.0000000
-7.0000000  -5.0000000  -2.0000000  -7.0000000
```

There are many linear algebra functions in Darwin, to obtain further details look at some of these help files: Cholesky CovarianceAnalysis CreateArray Eigenvalues GaussElim GivensElim Identity LinearRegression list matrix SvdAnalysis SvdBestBasis transpose

© 2005 by Gaston Gonnet, Informatik, ETH Zurich

Index of bio-recipes

Last updated on Fri Sep 2 14:45:47 2005 by GhG

!!! Dieses Dokument stammt aus dem ETH Web-Archiv und wird nicht mehr gepflegt !!!
!!! This document is stored in the ETH Web archive and is no longer maintained !!!