SVD

Opis

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:

  1. rešavanje linearnih jednačina,
  2. za kompresiju podataka,
  3. u algoritmu za preporučivanje,
  4. za konstrukciju frekvencije reči u dokumentima,
  5. za računanje $2-norme$ itd.

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:

  1. Kolone matrice $U$ su sopstveni vektori $AA^T$ - levi singularni vektori matrice $A$
  2. Kolone matrice $V$ su sopstveni vektori $A^TA$ - desni singularni vektori matrice $A$
  3. Dijagonalni elementi $\Sigma$ - singularne vrednosti matrice $A$
  4. Rang matrice A je broj $\sigma _i<>0$

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$

Unapređenja algoritma

Složenost izračunavanja

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.

Predlaganje unapređenja

Ovim naučnim radom se predlažu sledeće modifikacije i unapređenja:

  1. modifikacije sabiranja,
  2. rank-1 modifikacije,
  3. smanjenje kompleksnosti

Moguća redukcija

Ukoliko pretpostavimo da je redosled vrednosti $\sigma_{i}$ u neopadajućem poretku možemo eliminisati male vrednosti tako da ogranicimo $i < r < q$.
Tada je $U$ dimenzija $p \times r$, $\sigma$ vektor dimenzije $r$ i $V$ dimenzija $r \times r$ i ovim smo dobili tanku SVD ranka-r.

Modifikacija sabiranja

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$ Kada postavimo $U'^TKV'=S'$ tada $K$ proširuje prodprostore a $U'$ i $V'$ rotiraju $ $$\begin{bmatrix}U & P\end{bmatrix}$$ $ i $ $$\begin{bmatrix}V & Q\end{bmatrix}$$ $

Rank-1 modifikacije

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 računanje $K$ se svodi u slučaju proširivanja na što može biti izvršeno u $O(r^2)$, a u slučaju smanjivanja na se može rešiti bez $P$ i $Q=(b-Vn)/\sqrt{1-n^Tn}$ korišćenjem samo $i$-tog reda $V$


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


Vreme potrebno za racuanje SVD od X1 - 73.94300150871277
vreme potrebno za racuanje SVD od X2 - 73.95400047302246
vreme potrebno za racunanje inkrementalnog SVD 33.55100870132446

Diskusija o algoritmu

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.

Smanjivanje kompleksnosti

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 i zatim se u svakom koraku umesto da se rotiraju $U$ i $V$ rotiraju značajno manje $U'$ i $V'$

Reference:

  1. Mladen Nikolić, Anđelka Zečević, Naučno izračunavanje, 2017, http://ni.matf.bg.ac.rs/materijali/ni.pdf

  2. Matthew Brand, Fast low-rank modifications of the thin singular value decomposition, 2006, http://www.sciencedirect.com/science/article/pii/S0024379505003812