Mercurial > hg > segmentation

annotate pymf/svd.py @ 19:890cfe424f4a tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
added annotations
author mitian
date Fri, 11 Dec 2015 09:47:40 +0000
parents 26838b1f560f
children
rev   line source
mi@0 1 #!/usr/bin/python
mi@0 2 #
mi@0 3 # Copyright (C) Christian Thurau, 2010.
mi@0 4 # Licensed under the GNU General Public License (GPL).
mi@0 5 # http://www.gnu.org/licenses/gpl.txt
mi@0 6 """
mi@0 7 PyMF Singular Value Decomposition.
mi@0 8
mi@0 9 SVD : Class for Singular Value Decomposition
mi@0 10 pinv() : Compute the pseudoinverse of a Matrix
mi@0 11
mi@0 12 """
mi@0 13
mi@0 14
mi@0 15
mi@0 16 from numpy.linalg import eigh
mi@0 17 import scipy.sparse
mi@0 18
mi@0 19 try:
mi@0 20 import scipy.sparse.linalg.eigen.arpack as linalg
mi@0 21 except (ImportError, AttributeError):
mi@0 22 import scipy.sparse.linalg as linalg
mi@0 23
mi@0 24
mi@0 25 import numpy as np
mi@0 26
mi@0 27 def pinv(A, k=-1, eps=10**-8):
mi@0 28 # Compute Pseudoinverse of a matrix
mi@0 29 # calculate SVD
mi@0 30 svd_mdl = SVD(A, k=k)
mi@0 31 svd_mdl.factorize()
mi@0 32
mi@0 33 S = svd_mdl.S
mi@0 34 Sdiag = S.diagonal()
mi@0 35 Sdiag = np.where(Sdiag >eps, 1.0/Sdiag, 0.0)
mi@0 36
mi@0 37 for i in range(S.shape[0]):
mi@0 38 S[i,i] = Sdiag[i]
mi@0 39
mi@0 40 if scipy.sparse.issparse(A):
mi@0 41 A_p = svd_mdl.V.T * (S * svd_mdl.U.T)
mi@0 42 else:
mi@0 43 A_p = np.dot(svd_mdl.V.T, np.core.multiply(np.diag(S)[:,np.newaxis], svd_mdl.U.T))
mi@0 44
mi@0 45 return A_p
mi@0 46
mi@0 47
mi@0 48 class SVD():
mi@0 49 """
mi@0 50 SVD(data, show_progress=False)
mi@0 51
mi@0 52
mi@0 53 Singular Value Decomposition. Factorize a data matrix into three matrices s.t.
mi@0 54 F = | data - USV| is minimal. U and V correspond to eigenvectors of the matrices
mi@0 55 data*data.T and data.T*data.
mi@0 56
mi@0 57 Parameters
mi@0 58 ----------
mi@0 59 data : array_like [data_dimension x num_samples]
mi@0 60 the input data
mi@0 61
mi@0 62 Attributes
mi@0 63 ----------
mi@0 64 U,S,V : submatrices s.t. data = USV
mi@0 65
mi@0 66 Example
mi@0 67 -------
mi@0 68 >>> import numpy as np
mi@0 69 >>> data = np.array([[1.0, 0.0, 2.0], [0.0, 1.0, 1.0]])
mi@0 70 >>> svd_mdl = SVD(data, show_progress=False)
mi@0 71 >>> svd_mdl.factorize()
mi@0 72 """
mi@0 73
mi@0 74 _EPS=10**-8
mi@0 75
mi@0 76 def __init__(self, data, k=-1, rrank=0, crank=0):
mi@0 77 self.data = data
mi@0 78 (self._rows, self._cols) = self.data.shape
mi@0 79 if rrank > 0:
mi@0 80 self._rrank = rrank
mi@0 81 else:
mi@0 82 self._rrank = self._rows
mi@0 83
mi@0 84 if crank > 0:
mi@0 85 self._crank = crank
mi@0 86 else:
mi@0 87 self._crank = self._cols
mi@0 88
mi@0 89 # set the rank to either rrank or crank
mi@0 90 self._k = k
mi@0 91
mi@0 92 def frobenius_norm(self):
mi@0 93 """ Frobenius norm (||data - USV||) for a data matrix and a low rank
mi@0 94 approximation given by SVH using rank k for U and V
mi@0 95
mi@0 96 Returns:
mi@0 97 frobenius norm: F = ||data - USV||
mi@0 98 """
mi@0 99 if scipy.sparse.issparse(self.data):
mi@0 100 err = self.data - self.U*self.S*self.V
mi@0 101 err = err.multiply(err)
mi@0 102 err = np.sqrt(err.sum())
mi@0 103 else:
mi@0 104 err = self.data[:,:] - np.dot(np.dot(self.U, self.S), self.V)
mi@0 105 err = np.sqrt(np.sum(err**2))
mi@0 106
mi@0 107 return err
mi@0 108
mi@0 109
mi@0 110 def factorize(self):
mi@0 111 def _right_svd():
mi@0 112 AA = np.dot(self.data[:,:], self.data[:,:].T)
mi@0 113 values, u_vectors = eigh(AA)
mi@0 114
mi@0 115 # get rid of too low eigenvalues
mi@0 116 u_vectors = u_vectors[:, values > self._EPS]
mi@0 117 values = values[values > self._EPS]
mi@0 118
mi@0 119 # sort eigenvectors according to largest value
mi@0 120 idx = np.argsort(values)
mi@0 121 values = values[idx[::-1]]
mi@0 122
mi@0 123 # argsort sorts in ascending order -> access is backwards
mi@0 124 self.U = u_vectors[:,idx[::-1]]
mi@0 125
mi@0 126 # compute S
mi@0 127 self.S = np.diag(np.sqrt(values))
mi@0 128
mi@0 129 # and the inverse of it
mi@0 130 S_inv = np.diag(np.sqrt(values)**-1)
mi@0 131
mi@0 132 # compute V from it
mi@0 133 self.V = np.dot(S_inv, np.dot(self.U[:,:].T, self.data[:,:]))
mi@0 134
mi@0 135
mi@0 136 def _left_svd():
mi@0 137 AA = np.dot(self.data[:,:].T, self.data[:,:])
mi@0 138 values, v_vectors = eigh(AA)
mi@0 139
mi@0 140 # get rid of too low eigenvalues
mi@0 141 v_vectors = v_vectors[:, values > self._EPS]
mi@0 142 values = values[values > self._EPS]
mi@0 143
mi@0 144 # sort eigenvectors according to largest value
mi@0 145 # argsort sorts in ascending order -> access is backwards
mi@0 146 idx = np.argsort(values)[::-1]
mi@0 147 values = values[idx]
mi@0 148
mi@0 149 # compute S
mi@0 150 self.S= np.diag(np.sqrt(values))
mi@0 151
mi@0 152 # and the inverse of it
mi@0 153 S_inv = np.diag(1.0/np.sqrt(values))
mi@0 154
mi@0 155 Vtmp = v_vectors[:,idx]
mi@0 156
mi@0 157 self.U = np.dot(np.dot(self.data[:,:], Vtmp), S_inv)
mi@0 158 self.V = Vtmp.T
mi@0 159
mi@0 160 def _sparse_right_svd():
mi@0 161 ## for some reasons arpack does not allow computation of rank(A) eigenvectors (??) #
mi@0 162 AA = self.data*self.data.transpose()
mi@0 163 if self.data.shape[0] > 1:
mi@0 164 # do not compute full rank if desired
mi@0 165 if self._k > 0 and self._k < self.data.shape[0]-1:
mi@0 166 k = self._k
mi@0 167 else:
mi@0 168 k = self.data.shape[0]-1
mi@0 169
mi@0 170 values, u_vectors = linalg.eigen_symmetric(AA,k=k)
mi@0 171 else:
mi@0 172 values, u_vectors = eigh(AA.todense())
mi@0 173
mi@0 174 # get rid of too low eigenvalues
mi@0 175 u_vectors = u_vectors[:, values > self._EPS]
mi@0 176 values = values[values > self._EPS]
mi@0 177
mi@0 178 # sort eigenvectors according to largest value
mi@0 179 idx = np.argsort(values)
mi@0 180 values = values[idx[::-1]]
mi@0 181
mi@0 182 # argsort sorts in ascending order -> access is backwards
mi@0 183 self.U = scipy.sparse.csc_matrix(u_vectors[:,idx[::-1]])
mi@0 184
mi@0 185 # compute S
mi@0 186 self.S = scipy.sparse.csc_matrix(np.diag(np.sqrt(values)))
mi@0 187
mi@0 188 # and the inverse of it
mi@0 189 S_inv = scipy.sparse.csc_matrix(np.diag(1.0/np.sqrt(values)))
mi@0 190
mi@0 191 # compute V from it
mi@0 192 self.V = self.U.transpose() * self.data
mi@0 193 self.V = S_inv * self.V
mi@0 194
mi@0 195 def _sparse_left_svd():
mi@0 196 # for some reasons arpack does not allow computation of rank(A) eigenvectors (??)
mi@0 197 AA = self.data.transpose()*self.data
mi@0 198
mi@0 199 if self.data.shape[1] > 1:
mi@0 200 # do not compute full rank if desired
mi@0 201 if self._k > 0 and self._k < self.data.shape[1]-1:
mi@0 202 k = self._k
mi@0 203 else:
mi@0 204 k = self.data.shape[1]-1
mi@0 205 values, v_vectors = linalg.eigen_symmetric(AA,k=k)
mi@0 206 else:
mi@0 207 values, v_vectors = eigh(AA.todense())
mi@0 208 # get rid of too low eigenvalues
mi@0 209 v_vectors = v_vectors[:, values > self._EPS]
mi@0 210 values = values[values > self._EPS]
mi@0 211
mi@0 212 # sort eigenvectors according to largest value
mi@0 213 idx = np.argsort(values)
mi@0 214 values = values[idx[::-1]]
mi@0 215
mi@0 216 # argsort sorts in ascending order -> access is backwards
mi@0 217 self.V = scipy.sparse.csc_matrix(v_vectors[:,idx[::-1]])
mi@0 218
mi@0 219 # compute S
mi@0 220 self.S = scipy.sparse.csc_matrix(np.diag(np.sqrt(values)))
mi@0 221
mi@0 222 # and the inverse of it
mi@0 223 S_inv = scipy.sparse.csc_matrix(np.diag(1.0/np.sqrt(values)))
mi@0 224
mi@0 225 self.U = self.data * self.V * S_inv
mi@0 226 self.V = self.V.transpose()
mi@0 227
mi@0 228
mi@0 229 if self._rows > self._cols:
mi@0 230 if scipy.sparse.issparse(self.data):
mi@0 231 _sparse_left_svd()
mi@0 232 else:
mi@0 233 _left_svd()
mi@0 234 else:
mi@0 235 if scipy.sparse.issparse(self.data):
mi@0 236 _sparse_right_svd()
mi@0 237 else:
mi@0 238 _right_svd()
mi@0 239
mi@0 240 if __name__ == "__main__":
mi@0 241 import doctest
mi@0 242 doctest.testmod()