fml  0.1-0
Fused Matrix Library
stats.hh
1 // This file is part of fml which is released under the Boost Software
2 // License, Version 1.0. See accompanying file LICENSE or copy at
3 // https://www.boost.org/LICENSE_1_0.txt
4 
5 #ifndef FML_CPU_FUTURE_STATS_H
6 #define FML_CPU_FUTURE_STATS_H
7 #pragma once
8 
9 
10 #include "../linalg.hh"
11 
12 
13 namespace stats
14 {
15  template <typename T>
16  T sgn(T x)
17  {
18  return (T>0?1:(T<0?-1:0));
19  }
20 
21 
22 
23  namespace
24  {
25  template <typename REAL>
26  void ambpc(cpumat<REAL> &a, const cpumat<REAL> &b, const cpumat<REAL> &c)
27  {
28  const len_t m = a.nrows();
29  const len_t n = a.ncols();
30 
31  REAL *a_d = a.data_ptr();
32  REAL *b_d = b.data_ptr();
33  REAL *c_d = c.data_ptr();
34 
35  for (len_t j=0; j<n; j++)
36  {
37  for (len_t i=0; i<m; i++)
38  a_d[i + m*j] = a_d[i + m*j] - b[i + m*j] + c[i + m*j];
39  }
40  }
41 
42  template <typename REAL>
43  void shrink_op(cpuvec<REAL> &x, const REAL tau)
44  {
45  const len_t n = x.size();
46  const REAL *x_d = x.data_ptr();
47 
48  for (len_t i=0; i<n; i++)
49  x_d[i] = sgn(x[i]) * std::max(fabs(x[i]) - tau, 0);
50  }
51 
52  template <typename REAL>
53  void shrink_op(cpumat<REAL> &x, const REAL tau)
54  {
55  const len_t m = x.nrows();
56  const len_t n = x.ncols();
57  const REAL *x_d = x.data_ptr();
58 
59  for (len_t j=0; j<n; j++)
60  {
61  for (len_t i=0; i<m; i++)
62  x_d[i + m*j] = sgn(x[i + m*j]) * std::max(fabs(x[i + m*j]) - tau, 0);
63  }
64  }
65 
66  template <typename REAL>
67  void sv_thresh(cpumat<REAL> &x, const REAL tau, cpuvec<REAL> &s, cpumat<REAL> &u, cpumat<REAL> &vt)
68  {
69  linalg::svd(x, s, u, vt);
70  shrink_op(s, tau);
71 
72  colprod(vt, s);
73  matmult(false, false, (REAL)1.0, u, vt, x);
74  }
75  }
76 
77  template <typename REAL>
78  void robpca(cpumat<REAL> &M, const REAL delta=1e-7, const uint32_t maxiter=1000)
79  {
80  const len_t n1 = M.nrows();
81  const len_t n2 = M.ncols();
82 
83  const REAL lambda = 1/sqrt((REAL)std::max(n1, n2))
84  const REAL mu = 0.25 * ((REAL)n1*n2) / norm_1(M);
85 
86  cpumat<REAL> S(n1, n2);
87  cpumat<REAL> Y(n1, n2);
88  cpumat<REAL> L(n1, n2);
89 
90  S.fill_zero();
91  Y.fill_zero();
92 
93  cpuvec<REAL> s;
94  cpumat<REAL> u;
95  cpumat<REAL> vt;
96 
97  bool conv = false;
98  uint32_t iter = 0;
99  REAL term;
100 
101  REAL ub = delta * norm_F(M);
102 
103  while (!conv && iter<maxiter)
104  {
105  copy::cpu2cpu(M, L);
106  if (iter > 0)
107  ambpc(L, S, Y);
108 
109  sv_thresh(L, 1/mu, s, u, vt);
110 
111  // S = M - L + Y
112  copy::cpu2cpu(M, S);
113  ambpc(S, L, Y);
114 
115  sv_thresh(S, lambda/mu, s, u, vt);
116 
117  // Y = M - L + Y
118 
119 
120 
121  term = norm_F(...);
122  conv = (term <= ub);
123  iter++;
124  }
125  }
126 }
127 
128 
129 #endif
fml::linalg::norm_1
REAL norm_1(const cpumat< REAL > &x)
Computes the 1 matrix norm, the maximum absolute column sum.
Definition: norm.hh:38
fml::linalg::matmult
void matmult(const bool transx, const bool transy, const REAL alpha, const cpumat< REAL > &x, const cpumat< REAL > &y, cpumat< REAL > &ret)
Computes ret = alpha*op(x)*op(y) where op(A) is A or A^T.
Definition: matmult.hh:43
fml::linalg::norm_F
REAL norm_F(const cpumat< REAL > &x)
Computes the Frobenius/Euclidean matrix norm.
Definition: norm.hh:118