1. PLU rastav¶

Implementacija funkcije koja računa $PLU$ rastav matrice A

In [6]:
import numpy as np
def PLU_dekompozicija(A):
    n=len(A)
    p= np.arange(n)
    for k in range(0,n-1):
        #pivotiranje
        m=k+np.argmax(A[k:n,k])
        p[[k,m]]=p[[m,k]]
        A[[k,m]]=A[[m,k]]
        for i in range(k+1,n):
            if A[i,k]!=0.0:
                A[i,k]=A[i,k]/(A[k,k]) #multiplikator 
                A[i, k+1:n]=A[i, k+1:n]-A[i,k]*A[k, k+1:n]
    return [A,p]

Funkciju ćemo primijeniti na matrici $$\begin{bmatrix} 3& 17& 10\\ 2& 4& -2\\ 6& 18& -12\end{bmatrix}$$ Kao izlaz dobivamo matricu čiji gornji trokut je matrica $U$ a elementi ispod glavne dijagonale su elementi matrice $L$

In [7]:
A=np.array([[3.0, 17.0, 10.0 ],[2.0, 4.0, -2.0], [6.0, 18.0, -12.0]])
A, p =PLU_dekompozicija(A)
print("A=",A)
print("p=",p)

A= [[  6.          18.         -12.        ]
 [  0.5          8.          16.        ]
 [  0.33333333  -0.25         6.        ]]
p= [2 0 1]

Matricu $L$ možemo dobiti preko funkcije np.tril (proučiti kako ova funkcija radi), a matricu $U$ pomoću np.triu(ulazna matrica, id dijagonale)

Funkcija tril(ulazna matrica, identifikator dijagonale ) vraća matricu koja uzima donji dio ulazne matrice ispod(i uključujući) dijagonalu s time da glavna dijagonala ima identfikator 1, dijagonala iznad nje 1 i dijagonala ispod -1 itd. Slično je sa funkcijom triu(ulazna matrica, id dijagonale)

In [8]:
L=np.tril(A,-1)+np.identity(3)
print(L)

[[ 1.          0.          0.        ]
 [ 0.5         1.          0.        ]
 [ 0.33333333 -0.25        1.        ]]
In [9]:
U=np.triu(A,0)
print(U)

[[  6.  18. -12.]
 [  0.   8.  16.]
 [  0.   0.   6.]]
In [10]:
P_matrix=np.identity(3)[p]
print(P_matrix)

[[0. 0. 1.]
 [1. 0. 0.]
 [0. 1. 0.]]

Rješavanje sustava Ly=b (supstitucije unaprijed). u donjem kodu funkcija np.dot(x,y) predstavja skalarni produkt vektora. (Dijagonalni elementi od $L$ su 1 pa nije ni potrebni dijeliti s dijagonalnim elementima u donjem kodu)

In [11]:
def forward(L,b):
# supstitucija unaprijed
    b[0] = b[0] / L[0, 0]
    n=len(L)
    for i in range(1, n):
        b[i] = (b[i] - np.dot(L[i,:i], b[:i])) / L[i,i]
    return b
In [12]:
b=np.array([[67],[4],[6]])
In [13]:
y=forward(L,b[p])
print(y)

[[ 6]
 [64]
 [18]]

Rješavanje sustava $Ux=y$ odnosno povratne supstitucije

In [14]:
def backward(U,b):
# supstitucija unatrag
    n=len(U)
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - np.dot(U[k,k+1:n],b[k+1:n]))/U[k,k]
    return b
In [15]:
x=backward(U,y)
print(x)

[[1]
 [2]
 [3]]

Dobivanje PLU faktorizacije pomoću Pythonove biblioteke SciPy¶

lu funkcija iz scipy.linalg paketa kao izlaz daje matrice $P$, $L$, $U$, tako da je $A=PLU$, dakle matrica $P$ koju daje ova funkcija je zapravo transponirana matrica $P$ koju smo uveli na predavanjima, jer nam je dobro poznato $P^TP=I$. Dakle imamo: $P^T A=LU$

In [20]:
from scipy.linalg import lu
In [21]:
A=np.array([[3.0, 17.0, 10.0 ],[2.0, 4.0, -2.0], [6.0, 18.0, -12.0]])
In [22]:
P, L, U = lu(A)
In [23]:
print("P=\n", P)
print("L=\n", L)
print("U=\n", U)

P=
 [[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]
L=
 [[ 1.          0.          0.        ]
 [ 0.5         1.          0.        ]
 [ 0.33333333 -0.25        1.        ]]
U=
 [[  6.  18. -12.]
 [  0.   8.  16.]
 [  0.   0.   6.]]
In [24]:
print("Permutacijska matrica=\n",np.transpose(P))

Permutacijska matrica=
 [[0. 0. 1.]
 [1. 0. 0.]
 [0. 1. 0.]]

Faktorizacija Choleskog (za pozitivno definitne matrice)¶

In [25]:
import numpy as np
from numpy import dot
from math import sqrt

def cholesky(A):
    n = len(A)
    for k in range(n):
        A[k,k] = sqrt(A[k,k] - dot(A[k,0:k],A[k,0:k]))
        for i in range(k+1,n):
            A[i,k] = (A[i,k] - dot(A[i,0:k],A[k,0:k]))/A[k,k]
    for k in range(1,n): A[0:k,k] = 0.0
    return A

Funkciju ćemo primijeniti na matrici $$\begin{bmatrix} 1& 2& 1\\ 2& 5& 2\\ 1& 2& 10\end{bmatrix}$$ Kao izlaz dobivamo matricu čiji gornji trokut je matrica $U$ a elementi ispod glavne dijagonale su elementi matrice $L$

In [28]:
A=np.array([[1.0, 2.0, 1.0 ],[2.0, 5.0, 2.0], [1.0, 2.0, 10.0]])
cholesky(A)
print("A=",A)

A= [[1. 0. 0.]
 [2. 1. 0.]
 [1. 0. 3.]]

Možemo i upotrijebiti funkciju definiranu unutar numpy biblioteke za dobivanje faktorizacije Choleskog.

In [31]:
A=np.array([[1.0, 2.0, 1.0 ],[2.0, 5.0, 2.0], [1.0, 2.0, 10.0]])
L = np.linalg.cholesky(A)
L
Out[31]:
array([[1., 0., 0.],
       [2., 1., 0.],
       [1., 0., 3.]])

Možemo vidjeti što su svojstvene vrijednosti ove matrice

In [32]:
eig=np.linalg.eigvalsh(A)
print("Eigenvalues=", eig)

Eigenvalues= [ 0.16876617  4.86131375 10.96992008]

Vidimo da smo dobili pozitivne svojstvene vrijednosti, te znamo da je matrica pozitivno definitna i mogli smo izračunati faktorizaciju Choleskog? Što se dogodi kada matrica nije pozitivno definitna? Ispišite onda njezine svojstvene vrijednosti.