Jedna od najkorisnijih dekompozicija je dekompozicija singularnih vrednosti (SVD - singular value decomposition). Ona nam pruža mogućnost rešavanja ili razumevanja problema koji su definisani sistemima jednačina zasnovanih na sigularnim ili blisko singularnim matricama. Dok druge metode ne pokazuju dobre rezultate u tim slučajevima.
Ovaj metod ima široku primenu:
Ne zahteva da matrica $A$ koju faktorišemo bude kvadratna.
Za $A$ dimenzija $p \times q$ postoji ortogonalna matrica $U$ dimenzija $p \times p$, dijagonalna matrica $\Sigma$ dimenzija $p \times q$ i ortogonalna matrica $V$ dimenzija $q \times q$
$$A=U \Sigma V^T$$
Predstavlja formulu SVD, međutim moguće je predstaviti je i sledećom formulom:
$$X=U' \Sigma' V^T$$
gde su $U$ dimenzija $p \times q$, $\Sigma$ dimenzija $q \times q$ i $V$ dimenzija $q \times q$. Ovakva formula predstavlja tanku SVD.
Sledeća slika ilistruje formule.
$\Sigma$ se u računarstvu najčešće posmatra kao vektor dužine $q$, a formula može biti predstavljena kao suma matrica ranka-1 $\sum_{i}u_{i}\sigma_{i}v_{i}^T$, gde $\sigma_{i}$ predstavljaju singularne vrednosti (sa diagonale $\Sigma$ ), $u_{i}$ i $v_{i}$ su singularni vektori odnosno kolone $U$ i $V$
Važi da su:
Na osnovu prethodnog lako se može doći do sledećih formula:
$A^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T$
$AV=U\Sigma$
$det(A^TA-\lambda I)=0$ odavde je moguće izračunati $V$ i $\Sigma$, a zatim uvrštavanjem u prethodnu formulu i $U$
Izračunavanje potpune singularne dekompozicije je računarski izuzetno zahtevno i potrebno je $\Theta(pq^2)$ operacija. Cilj je komplseksnost smanjiti na $O(pqr)$. Jako brzo nakon pojavljivanja računara i praktičnog SVD algoritma počelo tragati za efikasnijim metodima. Tako je nastala tanka SVD modifikacija i slične metode koje su se odnosile na promenu podataka u matrici.
Ovim naučnim radom se predlažu sledeće modifikacije i unapređenja:
Ako je $X=USV^T$ onda za sabiranje želimo da postavimo modifikaciju tako da: $X+AB^T=$$\begin{bmatrix}U & A\end{bmatrix}$$ $$\begin{bmatrix} S & 0 \\ 0 & I \end{bmatrix}$$ $$\begin{bmatrix} V && B \end{bmatrix}$$^T$. U slučaju kad je $rank(X+AB^T) \leqslant r+c<min(p,q)$, matrice $U, V, A, B$ su visoke tanke.
Dalje ako je $P$ ortogonalna baza $(I-UU^T)A$ i $R_A=P^T(I-UU^T)A$ onda
$ $$\begin{bmatrix}U & A\end{bmatrix}$$ = $$\begin{bmatrix} U & P\end{bmatrix}$$ $$\begin{bmatrix}I & U^TA \\ 0 & R_A\end{bmatrix}$$ $
Slično je za $QR_B=(I-VV^T)B$.
Kombinacijom ovih formula dobijamo $X+AB^T=$$\begin{bmatrix}U & P\end{bmatrix}$$ K $$\begin{bmatrix}V & Q\end{bmatrix}$$^T$
Praktično se koristi tabela operacija sa SVD vrednostima koje znamo da izračunamo za tu operaciju, kao traženu koju možemo postići zamenom odgovarajućih $a$ i $b^T$
Za proširivanje: $a=c, b^T=[0,0,...,1]$
Za smanjivanje: $a=-c, b^T=[0,0,...,1]$
Za zamenu: $a=d-c, b^T=[0,0,...,1]$
Za ponovno centriranje:
$a=X*(I-(1/q)ll^T), b^T=l^T=[1,1,...,1]$
Ako uzmemo da važi
In [1]:
def svd_upd(V, c):
#prosirujemo V
#V = np.vstack([V, np.zeros(V.shape[1])])
#kreiramo b i punimo 0
b = np.zeros(V.shape[0])
#dodajemo 1 na kraj
b[-1] = 1
#transponujemo
b = np.reshape(b, (b.shape[0], 1))
#punimo a
a = np.reshape(c, (-1, 1))
return a,b
def svd_down(V, X):
#kreiramo b i punimo 0
b = np.zeros(V.shape[0])
b[-1] = 1
#transponujemo
b = np.reshape(b, (b.shape[0], 1))
#punimo a
a = np.reshape(np.multiply(X[:,-1], -1), (-1, 1))
return a,b
def svd_rev(V,X, c):
#prosirujemo V
V = np.vstack([V, np.zeros(V.shape[1])])
#kreiramo b i punimo 0
b = np.zeros(V.shape[0])
b[-1] = 1
#transponujemo
b = np.reshape(b, (b.shape[0], 1))
#punimo a
a = np.reshape(X[:,-1] - c, (-1, 1))
return a,b
def svd_recenter(V, X):
#kreiramo b i punimo 1
ones = np.ones(V.shape[1])
b = np.reshape(ones, (-1, 1))
#parametri potrebni za a
n = np.reshape(np.dot(np.transpose(V), b), (-1, 1))
q = b - np.dot(V, n)
#punimo a
a = np.reshape(np.multiply((-1/q), np.dot(X, b)), (-1, 1))
return a,b
In [2]:
def rediagonalization(U,S,V,a,b):
m = np.reshape(np.dot(np.transpose(U), a), (-1, 1))
p = np.reshape(a - np.dot(U, m), (-1, 1))
Ra = np.linalg.norm(p)
P = np.reshape(np.multiply((1 / Ra), p), (-1, 1))
n = np.reshape(np.dot(np.transpose(V), b), (-1, 1))
q = b - np.dot(V, n)
Rb = np.linalg.norm(q)
Q = np.reshape(np.multiply((1 / Rb), q), (-1, 1))
k = S
K = np.zeros((k.shape[0] + 1, k.shape[0] + 1))
K[:-1,:-1] = k
stack = np.vstack(np.append(m, Ra))
t = np.reshape(np.append(n, Rb), (1, -1))
dot = np.dot(stack, t)
K = np.add(K, dot)
return K
In [3]:
#op predstavlja operaciju koja se izvrsava
#0-upd,1-dwn,2-rev,3-rec
def svd(U,S,V,X,c=None,op=0):
if op==0:
if type(c)==type(np.array([])):
a, b = svd_upd(V, c)
else:
return None,None,None
elif op==1:
a, b = svd_down(V, X)
elif op==2:
if type(c)==type(np.array([])):
a, b = svd_rev(V, X, c)
else:
return None,None,None
elif op==3:
a, b = svd_recenter(V, X)
else:
return None,None,None
k=rediagonalization(U,S,V,a,b)
Sn, Vn=np.linalg.eig(k)
Sn=np.diag(Sn)
Un=np.transpose(np.linalg.inv(Vn))
return Un, Sn, Vn
In [14]:
import numpy as np
import time
import scipy.linalg as lin
size=30000
a = np.random.randint(100,size=(size, size//10))
X1 = a#lin.orth(a)
start1 = time.time()
U1, S1, V1 = np.linalg.svd(X1, full_matrices = False)
end1 = time.time()
print('Vreme potrebno za racuanje SVD od X1 - ' + str(end1 - start1))
c = np.random.randint(100,size=(size, 1))
X2 = np.hstack((X1,c))
start2 = time.time()
U2, S2, V2 = np.linalg.svd(X2, full_matrices = False)
end2 = time.time()
print('vreme potrebno za racuanje SVD od X2 - ' + str(end2 - start2))
#U2.dot(np.diag(S2)).dot(V2), X2
start3 = time.time()
Un, Sig, Vn = svd(U1,S1,V1,X1,c,op=0)
end3 = time.time()
print('vreme potrebno za racunanje inkrementalnog SVD ' + str(end3 - start3))
#Vn.T.dot(Sig).dot(Un), X2
Prikazana metrika - vreme potrebno za izračunavanje predstavlja značajnu uštedu. Međutim prveliki problem je u tačnosti izračunavanja. Ne dobija se matrica dekompozicija kojom se može rekonstruisati polazna matrica. Prikazani algoritam je direktno izveden iz formula iznetih u naučnom radu [1] međutim očigledno je potrebno izvršiti modifikaciju. Citiranost rada od u oko 300 naučnih radova čini ga osnovom inkrementalnog SVD, međutim analizom implementacija na koje sam nailazio (u drugim programskim jezicima) nisam uspeo doći do tačnih rezultata.
Korišćenjem odogvarajuće implementacije prethodnih koraka, moguće je kompleksnost dodatno smanjiti. Tako je za operaciju dodavanja moguće postići umesto $O(p(r+c)^2)$ - $O(pr)$, za rediagojalizaciju umesto $O((r+c)^3)$ na $O(r^2)$, a za rotiranje bodprostora sa $O((p+q)(c+r)^2)$ na $O(r^3)$.
Takva unapređenja se postižu umesto prethodne formule koristi pojednostavljena
Mladen Nikolić, Anđelka Zečević, Naučno izračunavanje, 2017, http://ni.matf.bg.ac.rs/materijali/ni.pdf
Matthew Brand, Fast low-rank modifications of the thin singular value decomposition, 2006, http://www.sciencedirect.com/science/article/pii/S0024379505003812