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[2];
[-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[1] = alpha[1]*Vt[1];
[0.2973, -0.3337, -0.2059] = [0.2973, -0.3337, -0.2059]
> Vt[2]*Vt[2];
1.0000
See Also:
?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
Last updated on Fri Sep 2 14:45:47 2005 by GhG