#------------------------------------#
#                                    #
#    Week2: Rechnen mit Matrizen     #
#                                    #
#------------------------------------#


Die wichtigsten Befehle zum Anlegen von Matrizen sind:
   
 - matrix()
 - rbind() und cbind()       # rowbind, columnbind
 - outer()
 - diag()

und die wichtigsten Operationen sind:

 - %*%        fuer Matrix-Multiplikation
 - solve()    fuer Matrix-Inversion
 - t()        Transponieren
 - det()      Determinante
 - eigen()    Eigenwerte und Eigenvektoren
   
Schauen wir uns das an:
   
   
### Start R-Session ###
 
x = 1:15                       # Vektor
x

A = matrix(x)
A                                         # 15 x 1 Matrix
A = matrix(x,nrow=3,ncol=5)
A                                         # 3 x 5 Matrix
A = matrix(x,nrow=3,ncol=5,byrow=TRUE)
A                                         # zeilenweises Auffuellen

# Matrix-Elemente:
A[1,1]                         # die Indizierung faengt 
A[1,2]                         # immer bei 1 an
A[2,5]
A[2,]                          # 2.Zeile
A[,4]                          # 4.Spalte
A[,3:5]                        # 3.-5.Spalte
A[c(1,3),c(2,4)]               # 1. und 3. Zeile und 2. und 4.Spalte                   


v1 = c(1,1,1,1)
v2 = c(1,2,3,4)
B = rbind(v1,v2)                          # zeilenweises Zusammensetzen
B
C = cbind(v1,v2)                          # spaltenweises Zusammensetzen
C
BB = cbind(B,B)
BB

a = 1:3
a
b = c(4,4,4)
b
outer(a,b)
outer(a,b,FUN="+")
outer(a,b,Fun="+")             # Fehlermeldung, Gross- und 
                               # Kleinschreibung ist wichtig
a = 1:3
b = 1:6
outer(a,b)
outer(a,b,FUN="+")
outer(a,b,FUN="*")             # ist die default-Einstellung

# n x n - Einheitsmatrix:
n = 10
e = rep(1,n)
e
E = diag(e)
E

# n x n - Nullmatrix:
N = matrix(0,n,n)
N
NN = diag(rep(0,n))           # geht auch
NN

A
diag(A)                       # diag(Matrix) = Vektor
diag(diag(A))                 # diag(Vektor) = Matrix


# Rechenoperationen mit Matrizen:
#
# Generelle Regel fuer das Rechnen 
# mit Vektoren und Matrizen in R:
#
# R versucht immer, ELEMENTWEISE zu rechnen:

x = 1:5
x
1/x 
x - 7
x^2
exp(x)

# fuer Matrizen gilt dasselbe:
A
1/A
A - 7
A^2
A*A                         # das ist keine Matrix-Multiplikation
exp(A)                      # das ist kein Matrix-Exponential, sondern 
                            # elementweise: exp(a11), exp(a12), exp(a13),...

# Damit stellt sich die Frage: Was ist die Syntax 
# fuer Matrix-Multiplikation und das Matrix-Inverse?

B
BT = t(B)
BT                          # B transponiert

# Matrix-Multiplikation: %*%

M1 = B %*% BT
M1
M2 = BT %*% B
M2

det(M1)                     # Determinante
det(M2)     

# Matrix-Inverse: solve()

M1inv = solve(M1)
M1inv
M1inv %*% M1                # Einheitsmatrix, up to numerical noise
M1 %*% M1inv 
1/M1                        # NICHT das Inverse
M1 * 1/M1
M1 %*% (1/M1) 
M1 %*% 1/M1                 # Fehlermeldung

det(M2)                     # = 0, nicht invertierbar
solve(M2)                   # Fehlermeldung


# Obere/Untere Dreiecks-Matrizen:
M2
lower.tri(M2)
ML = M2 * lower.tri(M2)     # Zahl * TRUE/FALSE kann man machen,
ML                          # TRUE->1, FALSE->0
MU = M2 * upper.tri(M2)
MU
d = diag(M2)
d
D = diag(d)
D
M = ML+MU+D
M
M == M2



# EIGENWERTE und EIGENVEKTOREN

# Wenn Sie etwas groessere Programme in R schreiben wollen, 
# wollen Sie natuerlich nicht Zeile fuer Zeile in das Consolen- 
# Fenster eingeben, wo ja, wenn Sie return gedrueckt haben, die 
# Zeile, der Code, dann verschwunden ist, Sie sehen das dann ja  
# nicht mehr. 
# Also Sie brauchen ein anderes Fenster, ein Programm-Editor, 
# wo Sie erstmal Ihren ganzen Code reinschreiben koennen und 
# dann ganz am Ende wollen Sie das dann executen. 
#
# Dazu, wenn Sie wie ich in Windows mit dem RGui arbeiten, 
# muessen Sie auf 
#                    File->New script
#  
# gehen, und dann oeffnet sich ein neues Fenster. Dort koennen 
# Sie Ihren R-Code reinschreiben, und wenn Sie den Code dann 
# ausfuehren wollen, selektieren Sie die entsprechenden Code-
# Zeilen mit der Maus (also so dass das dann blau hervorgehoben 
# oder blau unterlegt ist) und dann muessen Sie 
#
#                      die Taste F5 
#
# druecken. Dann springen genau die selektierten Code-Zeilen in 
# das R-Console-Fenster und werden dann also ausgefuehrt. 
#
# Oeffnen Sie jetzt ein neues Skript-Fenster, geben Sie den 
# folgenden Code ein und selektieren und executen Sie dann 
# mit der Taste F5:

n = 5                                # we repeat this for n=1000
x = rnorm(n*n)                       # n^2 Zufallszahlen
x
M = matrix(x,nrow=n,ncol=n)          # n x n Zufallsmatrix
M

# Eigenwerte und Eigenvektoren:
eigen(M)

# wir schreiben das Resultat in eine Variable:
res = eigen(M)

# Info zur Datenstruktur:
str(res)
names(res)

# das Resultat von eigen(M) ist eine Liste mit 2 Elementen, 
# res[[1]]=res$values ist vom Typ vector und res[[2]]=res$vectors 
# ist vom Typ Matrix.
# Listen sind im wesentlichen dazu da, um Objekte mit unterschied-
# lichen Datentypen in einem Objekt zusammenfassen zu koennen. 
# So ist der Rueckgabewert von R-Funktionen haeufig vom Typ Liste, 
# damit verschieden Arten von Information in einem Objekt zurueck-
# gegeben werden koennen. 

res[[1]]
res$values

res[[2]]
res$vectors


# Machen wir das jetzt nochmal fuer eine 1000 x 1000 Matrix 
# mit 1 Mio Zufallszahlen:

n = 1000
x = rnorm(n*n)
M = matrix(x,nrow=n,ncol=n)
res = eigen(M)

# das geht noch recht fix, messen wir mal die Zeit:
t0 = Sys.time()
res = eigen(M)
Sys.time()-t0                 # etwa 5 Sekunden auf meinem Tablet-PC

plot(res$values)

# Vermutung: die Eigenwerte sind gleichverteilt in 
# einem Kreis mit Radius sqrt(n).
# Checken wir das:

phi = seq( from=0 , to=2*pi , length=500 )
z = sqrt(n) * exp( 1i * phi )

plot(res$values)
lines( z , col="red" )             # fuegt Punkte zu einem bestehen Plot 
                                   # hinzu, im Linien-Format
points(z, type="l", col="green")   # ok, das geht auch

# probieren wir noch groessere n:
n = 2000
x = rnorm(n*n)
M = matrix(x,nrow=n,ncol=n)
t0 = Sys.time()
res = eigen(M)
Sys.time()-t0                 # etwa 32 Sekunden auf meinem Tablet-PC, 
                              # 2.18 sec auf meinem Geraet zu Hause


# wir wollen nur die Eigenwerte wissen, aber nicht die 
# Eigenvektoren:
t0 = Sys.time()
res = eigen(M, only.values=TRUE )
Sys.time()-t0                 # etwa 15 Sekunden auf meinem Tablet-PC, 
                              # 1.42 sec auf meinem Geraet zu Hause


# und nochmal das Bild:
plot(res$values)
phi = seq( from=0 , to=2*pi , length=500 )
z = sqrt(n) * exp( 1i * phi )
lines( z , col="red" )            # sehr huebsch..



# was gibt die eigen()-Funktion im Fall 
# von nicht-diagonalisierbaren Matrizen?

test = cbind(c(1,1),c(0,1))
test                              # nicht diagonalisierbar..
eigen(test)