#!/usr/bin/env python
# coding: utf-8

# # numpy.linalg: Lineare Algebra mit NumPy

# Viele grundlegende Methoden aus der linearen Algebra müssen wir nicht immer wieder mühsam selbst implementieren, sondern wir können auf Funktionen aus dem Paket numpy.linalg zurückgreifen. Hier ein paar nützliche Beispiele:

# In[1]:


import numpy as np
np.set_printoptions(legacy='1.25')


# In[2]:


A = np.array([[1,1,1],
              [1,2,2], 
              [1,2,3]])

b = np.array([[1], 
              [2], 
              [3]])


# ## Rang einer Matrix

# In[3]:


rk = np.linalg.matrix_rank(A)
rk


# ## Determinante einer Matrix

# In[4]:


detA = np.linalg.det(A)
detA


# ## Norm eines Vektors

# Berechnung der Euklidischen Norm eines Vektors $x\in\mathbb{C}^n$:
# \begin{align}
# \|x\| = \sqrt{\sum_{i=1}^{n} |x_i|^2} 
# \end{align}

# In[5]:


# per "Hand"
n1 = np.sqrt(sum(b*b))
n1[0]


# In[6]:


# mit numpy.linalg
n2 = np.linalg.norm(b)
n2


# ## Inverse einer Matrix

# In[7]:


Ainv = np.linalg.inv(A)
Ainv


# In[8]:


# Probe:
A @ Ainv


# In[9]:


np.linalg.norm(A @ Ainv - np.eye(3))


# ## Eigenwerte und Eigenvektoren

# In[10]:


ew, ev = np.linalg.eig(A)


# In[11]:


ew # i-ter Eigenwert an i-ter Stelle


# In[12]:


ev # Eigenvektor zum i-ten Eigenwert in i-ter Spalte


# Sind die Eigenvektoren schon normiert?

# In[13]:


[np.linalg.norm(ev[:,j]) for j in range(ev.shape[0])]


# Wird die Eigenwert-Gleichung $Av = \lambda v$ erfüllt?

# In[14]:


# Probe für den ersten Eigenvektor:
A @ ev[:,0] - ew[0] * ev[:,0]


# In[15]:


np.linalg.norm(A @ ev[:,0] - ew[0] * ev[:,0])


# In[16]:


# Probe für alle Eigenwerte und Eigenvektoren gleichzeitig:
np.linalg.norm(A @ ev - ev @ np.diag(ew))


# In[17]:


# Orthogonalität von ev
np.linalg.norm(ev @ ev.T - np.eye(3))


# In[18]:


# da alle Eigenwerte verschieden:
# Spektralzerlegung A = U D U^T
np.linalg.norm(A - ev @ np.diag(ew) @ ev.T)


# Wenn man nur die Eigenwerte benötigt, kann man auch numpy.linalg.eigvals() benutzen. Für symmetrisch bzw. hermitesche Matrizen gibt es die Funktion numpy.linalg.eigh().

# ## Lösung von $Ax = b$

# Für gegebenen Vektor b und invertierbare Matrix A passender Dimensionen wird ein Vektor x gesucht, sodass $Ax=b$ gilt.

# In[19]:


# Variante 1 ("Hau drauf"-Methode, NICHT empfohlen!)
x = np.linalg.inv(A) @ b
x


# In[20]:


# Variante 2 (empfohlen)
x = np.linalg.solve(A,b)
x


# In[21]:


# Probe:
np.linalg.norm(A@x - b)


# ## Sonstiges

# In[22]:


# Manche mögen es auch numpy.linalg einen eigenen Namen zu geben:
import numpy.linalg as la


# In[23]:


# Dann werden obige Befehle alle etwas kürzer, z.B.
la.det(A)


# In[24]:


# Alternative Überprüfung, ob z.B. Eigenwertgleichung erfüllt ist:
A @ ev[:,0] == ew[0] * ev[:,0] # hier nicht zufriedenstellend


# In[25]:


np.allclose(A @ ev[:,0], ew[0] * ev[:,0]) # hier gute Alternative


# In[26]:


# Probleme mit ganzzahligen Matrizen
A


# In[27]:


A[1,0] = 2.2


# In[28]:


A # Eintrag ist 1 statt 1.2 (es wurde also gerundet)


# In[29]:


# Variante 1
A = A.astype('float')
A


# In[30]:


A[1,0] = 2.2
A


# In[31]:


# besser:
A = np.array([[1,2,2,3,1], [2,4,4,6,2], [3,6,6,9,6], [1,2,4,5,3]], dtype='float')
A


# In[32]:


A[1,0] = 2.2
A


# In[ ]: