5 #ifndef FML_CPU_FUTURE_STATS_H
6 #define FML_CPU_FUTURE_STATS_H
10 #include "../linalg.hh"
18 return (T>0?1:(T<0?-1:0));
25 template <
typename REAL>
26 void ambpc(cpumat<REAL> &a,
const cpumat<REAL> &b,
const cpumat<REAL> &c)
28 const len_t m = a.nrows();
29 const len_t n = a.ncols();
31 REAL *a_d = a.data_ptr();
32 REAL *b_d = b.data_ptr();
33 REAL *c_d = c.data_ptr();
35 for (len_t j=0; j<n; j++)
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];
42 template <
typename REAL>
43 void shrink_op(cpuvec<REAL> &x,
const REAL tau)
45 const len_t n = x.size();
46 const REAL *x_d = x.data_ptr();
48 for (len_t i=0; i<n; i++)
49 x_d[i] = sgn(x[i]) * std::max(fabs(x[i]) - tau, 0);
52 template <
typename REAL>
53 void shrink_op(cpumat<REAL> &x,
const REAL tau)
55 const len_t m = x.nrows();
56 const len_t n = x.ncols();
57 const REAL *x_d = x.data_ptr();
59 for (len_t j=0; j<n; j++)
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);
66 template <
typename REAL>
67 void sv_thresh(cpumat<REAL> &x,
const REAL tau, cpuvec<REAL> &s, cpumat<REAL> &u, cpumat<REAL> &vt)
69 linalg::svd(x, s, u, vt);
73 matmult(
false,
false, (REAL)1.0, u, vt, x);
77 template <
typename REAL>
78 void robpca(cpumat<REAL> &M,
const REAL delta=1e-7,
const uint32_t maxiter=1000)
80 const len_t n1 = M.nrows();
81 const len_t n2 = M.ncols();
83 const REAL lambda = 1/sqrt((REAL)std::max(n1, n2))
84 const REAL mu = 0.25 * ((REAL)n1*n2) /
norm_1(M);
86 cpumat<REAL> S(n1, n2);
87 cpumat<REAL> Y(n1, n2);
88 cpumat<REAL> L(n1, n2);
101 REAL ub = delta *
norm_F(M);
103 while (!conv && iter<maxiter)
109 sv_thresh(L, 1/mu, s, u, vt);
115 sv_thresh(S, lambda/mu, s, u, vt);