1
2 """
3 The :mod:`parsimony.algorithms.nipals` module includes several algorithms
4 that minimises an implicit loss function based on the NIPALS algorithm.
5
6 Algorithms may not store states. I.e., if they are classes, do not keep
7 references to objects with state in the algorithm objects. It should be
8 possible to copy and share algorithms between e.g. estimators, and thus they
9 should not depend on any state.
10
11 Created on Thu Feb 20 17:46:17 2014
12
13 Copyright (c) 2013-2014, CEA/DSV/I2BM/Neurospin. All rights reserved.
14
15 @author: Tommy Löfstedt
16 @email: lofstedt.tommy@gmail.com
17 @license: BSD 3-clause.
18 """
19 import numpy as np
20
21 try:
22 from . import bases
23 except ValueError:
24 import parsimony.algorithms.bases as bases
25 import parsimony.utils.maths as maths
26 import parsimony.utils.consts as consts
27 import parsimony.utils.start_vectors as start_vectors
28 import parsimony.functions.penalties as penalties
29
30 __all__ = ["FastSVD", "FastSparseSVD", "FastSVDProduct", "PLSR"]
31
32
33
34
35 -class FastSVD(bases.ImplicitAlgorithm):
36
38 """A kernel SVD implementation.
39
40 Performs SVD of given matrix. This is always faster than np.linalg.svd.
41 Particularly, this is a lot faster than np.linalg.svd when M << N or
42 M >> N, for an M-by-N matrix.
43
44 Parameters
45 ----------
46 X : Numpy array. The matrix to decompose.
47
48 max_iter : Non-negative integer. Maximum allowed number of iterations.
49 Default is 100.
50
51 eps : Positive float. The tolerance used by the stopping criterion.
52
53 start_vector : BaseStartVector. A start vector generator. Default is
54 to use a random start vector.
55
56 Returns
57 -------
58 v : The right singular vector of X that corresponds to the largest
59 singular value.
60
61 Examples
62 --------
63 >>> import numpy as np
64 >>> from parsimony.algorithms.nipals import FastSVD
65 >>>
66 >>> np.random.seed(0)
67 >>> X = np.random.random((10, 10))
68 >>> fast_svd = FastSVD()
69 >>> fast_svd.run(X)
70 array([[-0.3522974 ],
71 [-0.35647707],
72 [-0.35190104],
73 [-0.34715338],
74 [-0.19594198],
75 [-0.24103104],
76 [-0.25578904],
77 [-0.29501092],
78 [-0.42311297],
79 [-0.27656382]])
80 >>>
81 >>> np.random.seed(0)
82 >>> X = np.random.random((100, 150))
83 >>> fast_svd = FastSVD()
84 >>> v = fast_svd.run(X)
85 >>> us = np.linalg.norm(np.dot(X, v))
86 >>> s = np.linalg.svd(X, full_matrices=False, compute_uv=False)
87 >>> abs(np.sum(us ** 2.0) - np.max(s) ** 2.0) < 5e-13
88 True
89 >>>
90 >>> np.random.seed(0)
91 >>> X = np.random.random((100, 50))
92 >>> fast_svd = FastSVD()
93 >>> v = fast_svd.run(X)
94 >>> us = np.linalg.norm(np.dot(X, v))
95 >>> s = np.linalg.svd(X, full_matrices=False, compute_uv=False)
96 >>> abs(np.sum(us ** 2.0) - np.max(s) ** 2.0)
97 4.5474735088646412e-13
98 """
99 if start_vector is None:
100 start_vector = start_vectors.RandomStartVector(normalise=True)
101 M, N = X.shape
102 if M < 80 and N < 80:
103 _, _, V = np.linalg.svd(X, full_matrices=True)
104 v = V[[0], :].T
105 elif M < N:
106 K = np.dot(X, X.T)
107 t = start_vector.get_vector(X.shape[0])
108 for it in xrange(max_iter):
109 t_ = t
110 t = np.dot(K, t_)
111 t *= 1.0 / np.sqrt(np.sum(t ** 2.0))
112
113 if maths.norm(t_ - t) / maths.norm(t) < eps:
114 break
115
116 v = np.dot(X.T, t)
117 v *= 1.0 / np.sqrt(np.sum(v ** 2.0))
118
119 else:
120 K = np.dot(X.T, X)
121 v = start_vector.get_vector(X.shape[1])
122 for it in xrange(max_iter):
123 v_ = v
124 v = np.dot(K, v_)
125 v *= 1.0 / np.sqrt(np.sum(v ** 2.0))
126
127 if maths.norm(v_ - v) / maths.norm(v) < eps:
128 break
129
130 return v
131
132
134
135 - def run(self, X, max_iter=100, start_vector=None):
136 """A kernel SVD implementation for sparse CSR matrices.
137
138 This is usually faster than np.linalg.svd when density < 20% and when
139 M << N or N << M (at least one order of magnitude). When M = N >= 10000
140 it is faster when the density < 1% and always faster regardless of
141 density when M = N < 10000.
142
143 These are ballpark estimates that may differ on your computer.
144
145 Parameters
146 ----------
147 X : Numpy array. The matrix to decompose.
148
149 max_iter : Integer. Maximum allowed number of iterations.
150
151 start_vector : BaseStartVector. A start vector generator. Default is
152 to use a random start vector.
153
154 Returns
155 -------
156 v : Numpy array. The right singular vector.
157
158 Example
159 -------
160 >>> import numpy as np
161 >>> from parsimony.algorithms.nipals import FastSparseSVD
162 >>> np.random.seed(0)
163 >>> X = np.random.random((10,10))
164 >>> fast_sparse_svd = FastSparseSVD()
165 >>> fast_sparse_svd.run(X)
166 array([[ 0.3522974 ],
167 [ 0.35647707],
168 [ 0.35190103],
169 [ 0.34715338],
170 [ 0.19594198],
171 [ 0.24103104],
172 [ 0.25578904],
173 [ 0.29501092],
174 [ 0.42311297],
175 [ 0.27656382]])
176 """
177 if start_vector is None:
178 start_vector = start_vectors.RandomStartVector(normalise=True)
179 M, N = X.shape
180 if M < N:
181 K = X.dot(X.T)
182 t = start_vector.get_vector(X.shape[0])
183 for it in xrange(max_iter):
184 t_ = t
185 t = K.dot(t_)
186 t *= 1.0 / np.sqrt(np.sum(t ** 2.0))
187
188 a = float(np.sqrt(np.sum((t_ - t) ** 2.0)))
189 if a < consts.TOLERANCE:
190 break
191
192 v = X.T.dot(t)
193 v *= 1.0 / np.sqrt(np.sum(v ** 2.0))
194
195 else:
196 K = X.T.dot(X)
197 v = start_vector.get_vector(X.shape[1])
198 for it in xrange(max_iter):
199 v_ = v
200 v = K.dot(v_)
201 v *= 1.0 / np.sqrt(np.sum(v ** 2.0))
202
203 a = float(np.sqrt(np.sum((v_ - v) ** 2.0)))
204 if a < consts.TOLERANCE:
205 break
206
207 return v
208
209
211
212 - def run(self, X, Y, start_vector=None,
213 eps=consts.TOLERANCE, max_iter=100, min_iter=1):
214 """A kernel SVD implementation of a product of two matrices, X and Y.
215 I.e. the SVD of np.dot(X, Y), but the SVD is computed without actually
216 computing the matrix product.
217
218 Performs SVD of a given matrix. This is always faster than
219 np.linalg.svd when extracting only one, or a few, vectors.
220
221 Parameters
222 ----------
223 X : Numpy array with shape (n, p). The first matrix of the product.
224
225 Y : Numpy array with shape (p, m). The second matrix of the product.
226
227 start_vector : Numpy array. The start vector.
228
229 eps : Float. Tolerance.
230
231 max_iter : Integer. Maximum number of iterations.
232
233 min_iter : Integer. Minimum number of iterations.
234
235 Returns
236 -------
237 v : Numpy array. The right singular vector of np.dot(X, Y) that
238 corresponds to the largest singular value of np.dot(X, Y).
239
240 Example
241 -------
242 >>> import numpy as np
243 >>> from parsimony.algorithms.nipals import FastSVDProduct
244 >>> np.random.seed(0)
245 >>> X = np.random.random((15,10))
246 >>> Y = np.random.random((10,5))
247 >>> fast_svd = FastSVDProduct()
248 >>> fast_svd.run(X, Y)
249 array([[ 0.47169804],
250 [ 0.38956366],
251 [ 0.41397845],
252 [ 0.52493576],
253 [ 0.42285389]])
254 """
255 M, N = X.shape
256
257 if start_vector is None:
258 start_vector = start_vectors.RandomStartVector(normalise=True)
259 v = start_vector.get_vector(Y.shape[1])
260
261 for it in xrange(1, max_iter + 1):
262 v_ = v
263 v = np.dot(X, np.dot(Y, v_))
264 v = np.dot(Y.T, np.dot(X.T, v))
265 v *= 1.0 / np.sqrt(np.sum(v ** 2.0))
266
267 if np.sqrt(np.sum((v_ - v) ** 2.0)) < eps \
268 and it >= min_iter:
269 break
270
271 return v
272
273
274 -class PLSR(bases.ImplicitAlgorithm,
275 bases.IterativeAlgorithm):
276 """A NIPALS implementation for PLS regresison.
277
278 Parameters
279 ----------
280 max_iter : Non-negative integer. Maximum allowed number of iterations.
281 Default is 200.
282
283 eps : Positive float. The tolerance used in the stopping criterion.
284
285 Examples
286 --------
287 >>> from parsimony.algorithms.nipals import PLSR
288 >>> import numpy as np
289 >>> np.random.seed(42)
290 >>>
291 >>> X = np.random.rand(10, 10)
292 >>> Y = np.random.rand(10, 5)
293 >>> w = np.random.rand(10, 1)
294 >>> c = np.random.rand(5, 1)
295 >>> plsr = PLSR()
296 >>> w, c = plsr.run([X, Y], [w, c])
297 >>> w
298 array([[ 0.34682103],
299 [ 0.32576718],
300 [ 0.28909788],
301 [ 0.40036279],
302 [ 0.32321038],
303 [ 0.39060766],
304 [ 0.22351433],
305 [ 0.28643062],
306 [ 0.29060872],
307 [ 0.23712672]])
308 >>> c
309 array([[ 0.29443832],
310 [ 0.35886751],
311 [ 0.33847141],
312 [ 0.23526002],
313 [ 0.35910191]])
314 >>> C, S, W = np.linalg.svd(np.dot(Y.T, X))
315 >>> w_ = W[0, :].reshape(10, 1)
316 >>> w_ = -w_ if w_[0, 0] < 0.0 else w_
317 >>> w = -w if w[0, 0] < 0.0 else w
318 >>> round(np.linalg.norm(w - w_), 15)
319 1.52884e-10
320 >>> abs(np.dot(np.dot(X, w).T,
321 ... np.dot(Y, c / np.linalg.norm(c)))[0, 0] - S[0]) < 5e-15
322 True
323 """
329
330 - def run(self, XY, wc=None):
331 """A NIPALS implementation for PLS regresison.
332
333 Parameters
334 ----------
335 XY : List of two numpy arrays. XY[0] is n-by-p and XY[1] is n-by-q. The
336 independent and dependent variables.
337
338 wc : List of numpy array. The start vectors.
339
340 Returns
341 -------
342 w : Numpy array, p-by-1. The weight vector of X.
343
344 c : Numpy array, q-by-1. The weight vector of Y.
345 """
346 X = XY[0]
347 Y = XY[1]
348
349 n, p = X.shape
350
351 if wc is not None:
352 w_new = wc[0]
353 else:
354 maxi = np.argmax(np.sum(Y ** 2.0, axis=0))
355 u = Y[:, [maxi]]
356 w_new = np.dot(X.T, u)
357 w_new *= 1.0 / maths.norm(w_new)
358
359 for i in range(self.max_iter):
360 w = w_new
361
362 c = np.dot(Y.T, np.dot(X, w))
363 w_new = np.dot(X.T, np.dot(Y, c))
364 normw = maths.norm(w_new)
365 if normw > 10.0 * consts.FLOAT_EPSILON:
366 w_new *= 1.0 / normw
367
368 if maths.norm(w_new - w) < maths.norm(w) * self.eps:
369 break
370
371 self.num_iter = i
372
373 t = np.dot(X, w)
374 tt = np.dot(t.T, t)[0, 0]
375 c = np.dot(Y.T, t)
376 if tt > consts.TOLERANCE:
377 c *= 1.0 / tt
378
379 return w_new, c
380
381
382 -class SparsePLSR(bases.ImplicitAlgorithm,
383 bases.IterativeAlgorithm):
384 """A NIPALS implementation for Sparse PLS regresison.
385
386 Parameters
387 ----------
388 l : List or tuple of two non-negative floats. The Lagrange multipliers, or
389 regularisation constants, for the X and Y blocks, respectively.
390
391 penalise_y : Bool. Whether or not to penalise the Y block as well.
392
393 max_iter : Non-negative integer. Maximum allowed number of iterations.
394 Default is 200.
395
396 eps : Positive float. The tolerance used in the stopping criterion.
397
398 Examples
399 --------
400 >>> from parsimony.algorithms.nipals import SparsePLSR
401 >>> import numpy as np
402 >>> np.random.seed(42)
403 >>>
404 >>> X = np.random.rand(10, 10)
405 >>> Y = np.random.rand(10, 5)
406 >>> w = np.random.rand(10, 1)
407 >>> c = np.random.rand(5, 1)
408 >>> plsr = SparsePLSR(l=[4.0, 5.0])
409 >>> w, c = plsr.run([X, Y], [w, c])
410 >>> w
411 array([[ 0.32012726],
412 [ 0.31873833],
413 [ 0.15539258],
414 [ 0.64271827],
415 [ 0.23337738],
416 [ 0.54819589],
417 [ 0. ],
418 [ 0.06088551],
419 [ 0. ],
420 [ 0. ]])
421 >>> c
422 array([[ 0.1463623 ],
423 [ 0.66483154],
424 [ 0.4666803 ],
425 [ 0. ],
426 [ 0.5646119 ]])
427 """
430
431 super(SparsePLSR, self).__init__(max_iter=max_iter, **kwargs)
432
433 self.eps = max(consts.TOLERANCE, float(eps))
434
435 self.l = [max(0.0, float(l[0])),
436 max(0.0, float(l[1]))]
437
438 self.penalise_y = bool(penalise_y)
439
440 - def run(self, XY, wc=None):
441 """A NIPALS implementation for sparse PLS regresison.
442
443 Parameters
444 ----------
445 XY : List of two numpy arrays. XY[0] is n-by-p and XY[1] is n-by-q. The
446 independent and dependent variables.
447
448 wc : List of numpy array. The start vectors.
449
450 Returns
451 -------
452 w : Numpy array, p-by-1. The weight vector of X.
453
454 c : Numpy array, q-by-1. The weight vector of Y.
455 """
456 X = XY[0]
457 Y = XY[1]
458
459 n, p = X.shape
460
461 l1_1 = penalties.L1(l=self.l[0])
462 l1_2 = penalties.L1(l=self.l[1])
463
464 if wc is not None:
465 w_new = wc[0]
466 else:
467 maxi = np.argmax(np.sum(Y ** 2.0, axis=0))
468 u = Y[:, [maxi]]
469 w_new = np.dot(X.T, u)
470 w_new *= 1.0 / maths.norm(w_new)
471
472 for i in range(self.max_iter):
473 w = w_new
474
475 c = np.dot(Y.T, np.dot(X, w))
476 if self.penalise_y:
477 c = l1_2.prox(c)
478 normc = maths.norm(c)
479 if normc > consts.TOLERANCE:
480 c *= 1.0 / normc
481
482 w_new = np.dot(X.T, np.dot(Y, c))
483 w_new = l1_1.prox(w_new)
484 normw = maths.norm(w_new)
485 if normw > consts.TOLERANCE:
486 w_new *= 1.0 / normw
487
488 if maths.norm(w_new - w) / maths.norm(w) < self.eps:
489 break
490
491 self.num_iter = i
492
493
494
495
496
497
498
499 return w_new, c
500
501
502 if __name__ == "__main__":
503 import doctest
504 doctest.testmod()
505