ea.h

 #ifndef M_EA_H
 #define M_EA_H
 
-#include "../LIBS/Eigen/Dense"
 #include "../CORE/config/config.h"
-#include "../CORE/eigen_types.h"
+#include "svd.h"
+#include "ridge_regression.h"
+#include "../CORE/core_types.h"
+#include "../CORE/lapack.h"
 
 class EA {
     private:
-        Config &config;
+        const Config &config;
+        const SVD &svd;
+        const t_type &T;
         int verbose;
     public:
-        EA(Config &c):
+        EA(const Config &c, const SVD &svd, const t_type &T):
             config(c),
+            svd(svd),
+            T(T),
             verbose(c.get<int>("VERBOSE"))
     {}
 
-        int run  (const mat &Phi,const vec &T,
-                const mat &S, const mat &U,const mat &V,
-                double &alpha, double &beta) {
-
+        int run(double &alpha, double &beta) {
             // Phi = USV^T
-            size_t N=Phi.rows();
+            size_t N=svd.shapeA().first;
             size_t maxsteps=40;
             const double EPS2 = 5e-12;
             int convergence_status = 1;
@@ -31,40 +34,41 @@ class EA {
             size_t counter=0;
             double lambda = alpha/beta;
 
+            double *s = svd.getS();
+
             // square of the S matrix contains all eigenvalues
             // of Phi^T*Phi = (VSU^T)(USV^T) = VS^2V^T = S^2
-            const vec eval = S.array().pow(2);
-
+            double *eval = new double[svd.sizeS()];
+            for (size_t i=0; i<svd.sizeS(); ++i)
+                eval[i] = s[i]*s[i];
 
             if (verbose) printf("%9s  %9s  %9s  %9s  %9s  %9s\n","del alpha", "del beta", "alpha", "beta", "lambda", "gamma");
             // D is a filter factors matrix
-            mat D(S.rows(),S.rows());
-            D.setZero();
+            double *d = new double[svd.sizeS()];
+
             int outprec=config.get<int>("OUTPREC");
 
+            double *t = new double[svd.shapeA().first];
+            double *x = new double[svd.shapeA().second];
             while (test1>EPS2 || test2>EPS2) {
 
-                for (long int i=0; i<S.rows(); ++i) {
-                    double s = S(i,0);
-                    if (s>std::numeric_limits<double>::min())
-                        D(i,i) = (s)/(s*s+lambda);
-                    else
-                        D(i,i) = 0.0;
-                }
                 // regularised least square problem (Tikhonov Regularization):
-                vec m_n = V*D*U.transpose()*T;
+                t_type m_n((size_t)svd.sizeS());
+                RidgeRegression::solve(svd,T,m_n,lambda);
 
                 double gamma = 0.0;
-                for (long int j=0; j<S.rows(); j++) {
-                    // exlude eigenvalues which are smaller than ~1e-10
-                    // for numerical stability b/c matrix is not regularized
-                    if (beta*eval(j) > std::numeric_limits<double>::min()) {
-                        gamma += (beta*eval(j))/(beta*eval(j)+alpha);
+                for (size_t j=0; j<svd.sizeS(); j++) {
+                    if (beta*eval[j] > std::numeric_limits<double>::min()) {
+                        gamma += (beta*eval[j])/(beta*eval[j]+alpha);
                     }
                 }
-                alpha = gamma / m_n.dot(m_n);
+                alpha = gamma / (m_n*m_n);
+
 
-                double temp = (T-Phi*m_n).dot(T-Phi*m_n);
+                t_type Tpred = svd.predict(m_n,t,x);
+                double temp = (T-Tpred)*(T-Tpred);
+
+                //double temp = (T-U*S.asDiagonal()*V.transpose()*m_n).dot(T-U*S.asDiagonal()*V.transpose()*m_n);
                 beta = (static_cast<double>(N)-gamma)/temp;
 
                 test1 = fabs(alpha_old-alpha)/alpha;
@@ -82,8 +86,45 @@ class EA {
             }
             if (convergence_status && verbose) std::cout << "Number of steps for convergence: " + to_string(counter,outprec) << std::endl;
             if (!convergence_status && verbose) std::cout << "Max number of steps reached: " + to_string(counter,outprec) << std::endl;
+
+            delete [] eval;
+            delete [] d;
+
+            delete[] t;
+            delete[] x;
+
             return convergence_status;
+        };
+        Matrix get_sigma(double alpha, double beta) {
+            double *vt = svd.getVT();
+            double *s = svd.getS();
+            int n = static_cast<int>(svd.shapeVT().first);
+
+            Matrix Sigma((size_t)n,(size_t)n);
+            Sigma.setZero();
+            double *sigma = &Sigma.data()[0];
+
+            // Build P_inv matrix
+            for (size_t i=0; i<svd.sizeS(); ++i) {
+                Sigma(i,i) = 1.0/(beta*s[i]*s[i]+alpha);
+            }
+
+            //Sigma = V*P_inv*V.transpose();
+            char transa = 'T';
+            char transb = 'N';
+            double alpha_ = 1.0;
+            double beta_ = 0.0;
+
+            double *work = new double[n*n];
+
+            dgemm_(&transb, &transb, &n, &n, &n, &alpha_, sigma, &n, vt, &n,
+                &beta_, work, &n);
+
+            dgemm_(&transa, &transb, &n, &n, &n, &alpha_, vt, &n, work, &n,
+                &beta_, sigma, &n);
 
+            delete[] work;
+            return Sigma;
         }
 };
 #endif

ekm.h

 #ifndef EKM_H
 #define EKM_H
 
-#include "../LIBS/Eigen/Dense"
-#include "../CORE/eigen_types.h"
 #include "../CORE/core_types.h"
 
 template <typename K>
@@ -10,17 +8,17 @@ class EKM {
     private:
         K kernel;
     public:
-        Eigen::MatrixXd KK;
+        Matrix KK;
         template <typename T>
             EKM(T &t):
                 kernel(t)
     {}
-        void configure(aeds_type2 &basis) {
-            // KK is temp matrix
-            KK = Eigen::MatrixXd(basis.size(), basis.size());
+        void configure(Matrix &basis) {
+            KK = Matrix(basis.cols(), basis.cols());
 
-            for (long int i=0; i<KK.rows(); i++) {
-                for (long int j=0; j<KK.cols(); j++) {
+            // Compute lower diagonal part only
+            for (size_t i=0; i<KK.cols(); i++) {
+                for (size_t j=i+1; j<KK.cols(); j++) {
                     double val = kernel(basis.col(i),basis.col(j));
                     KK(i,j) = val;
                 }
@@ -28,21 +26,53 @@ class EKM {
 
             // KK is a kernel matrix
             // we compute cholesky for KK^-1
-            // then KK = L^-1
-            KK = KK.inverse();
-            Eigen::LLT<Eigen::MatrixXd> llt(KK);
-            KK = llt.matrixL();
-            KK.transposeInPlace();
+            inverse(KK);
+            cholesky(KK,'L');
+
+        }
+        void inverse(Matrix &A) {
+            double *a = &A.data()[0];
+            int n = A.rows();   // A is square
+            int *ipiv = new int[n];
+            int lwork = n*n;
+            double *work = new double[lwork];
+            int info;
+
+            dgetrf_(&n,&n,a,&n,ipiv,&info);
+            dgetri_(&n,a,&n,ipiv,work,&lwork,&info);
+
+            delete[] ipiv;
+            delete[] work;
+        }
+        void cholesky(Matrix &A, char uplo) {
+            // matrix must be either U or L
+            int n = A.cols(); // A is square
+            double *a = &A.data()[0];
+            int info;
+            dpotrf_(&uplo, &n, a, &n, &info);
         }
 
         void project(phi_type &Phi) {
-            for (long int i=0; i<Phi.rows(); ++i) {
-                Phi.row(i)= KK*Phi.row(i).transpose();
-            }
+            char transa = 'N';
+            char transb = 'N';
+            int m = Phi.rows();
+            int n = Phi.cols();
+            double *a = &Phi.data()[0];
+            double *b = &KK.data()[0];
+            double *c = new double[m*n];
+            double alpha = 1.0;
+            double beta = 0.0;
+            dgemm_(&transa, &transb, &m, &n, &n, &alpha, a, &m, b, &n,
+                    &beta, c, &m);
+
+            Phi = Matrix(c,m,n);
+
+            delete[] c;
+
         }
-        aed_type2 project(aed_type2 &aed, aeds_type2 &basis) {
-            aed_type2 temp(basis.size());
-            for (long int b=0; b<basis.cols(); b++) {
+        aed_type2 project(aed_type2 &aed, Matrix &basis) {
+            aed_type2 temp(basis.cols());
+            for (size_t b=0; b<basis.cols(); b++) {
                 temp(b) = kernel(aed,basis[b]);
             }
             return  KK*temp;

ols.h

+#ifndef OLS_H
+#define OLS_H
+
+
+#include "../CORE/core_types.h"
+
+class OLS {
+    public:
+        template <typename M, typename V>
+            static void solve(M &A, V &B, V &weights) {
+                // A will be destroyed on exit
+                // B will contain solution
+                int m = A.rows();  //IN
+                int n = A.cols();  //IN
+                int nrhs = 1; // number of colums in B (targets T) //IN
+                double *a = &A.data()[0]; //IN
+                int lda = m;
+                double *b = &B.data()[0];    // IN/OUT
+                int ldb = std::max(m,n);
+                double *s = new double[std::min(m,n)];   // size: rows
+                double rcond = 1e-8;   // 1 use machine precision
+                int rank;
+                double *work;
+                int lwork=-1; // QUERRY
+                int liwork = std::max(
+                        0,int(std::log2(std::min(m,n)/(25+1))) + 1)
+                    + 11*std::min(m,n);
+                int smlsiz = 25;
+                int nlvl = std::max( 0, int( std::log2( std::min(m,n)/(smlsiz+1) ) )+1 );
+                int *iwork = new int[3*std::min(m,n)*nlvl + 11*std::min(m,n)];
+                int info;
+                double wkopt;
+                dgelsd_(&m,&n,&nrhs,a,&lda,b,&ldb,s,&rcond,&rank,
+                        &wkopt, &lwork, iwork, &info);
+                lwork = (int)wkopt;
+                work = new double[lwork];
+                dgelsd_(&m,&n,&nrhs,a,&lda,b,&ldb,s,&rcond,&rank,
+                        work, &lwork, iwork, &info);
+
+                for (size_t i=0;i<n;++i)
+                    std::cout << B(i)<< " ";
+                std::cout<< std::endl;
+                weights = B.copy(A.cols());
+
+                delete[] s;
+                delete[] work;
+                delete[] iwork;
+            }
+};
+#endif

ridge_regression.h

+#ifndef RIDGE_REGRESSION_H
+#define RIDGE_REGRESSION_H
+
+#include "../CORE/core_types.h"
+#include "../CORE/maths.h"
+#include "../CORE/lapack.h"
+#include "svd.h"
+
+class RidgeRegression {
+    public:
+        template <typename V, typename W>
+            static void solve(const SVD &svd, V /*&*/T, W &weights,
+                    const double lambda) {
+
+                double *u = svd.getU();
+                double *s = svd.getS();
+                double *vt = svd.getVT();
+
+                double *d = new double[svd.sizeS()];
+
+                for (size_t i=0; i<svd.sizeS(); ++i) {
+                    double sv = s[i];
+                    if (sv>std::numeric_limits<double>::min())
+                        d[i] = (sv)/(sv*sv+lambda);
+                    else
+                        d[i] = 0.0;
+                }
+
+                int m = svd.shapeA().first;
+                int n = svd.shapeA().second;
+                int ldu = svd.shapeU().first;
+                int ucol = svd.shapeU().second;
+                int lda = m;
+                int ldvt = svd.shapeVT().first;
+
+                // SOLVE: V D U^T T
+
+                //STEP1: U^T T
+                char trans = 'T';
+                double alpha_ = 1.0;
+                double *x = &T.data()[0];
+                double *y = new double[lda];
+                int incx = 1;
+                double beta_ = 0.0;
+                dgemv_(&trans, &ldu, &ucol, &alpha_, u, &lda, x, &incx,
+                        &beta_, y, &incx);
+
+                // STEP2 : D (U^T T)
+                for (int r=0; r<ldvt; ++r) {
+                    y[r]*=d[r];
+                }
+
+                // STEP3: V (D U^T T)
+                trans='T';
+                dgemv_(&trans, &ldvt, &n, &alpha_, vt, &ldvt, y, &incx,
+                        &beta_, x, &incx);
+
+                weights = t_type(x, (size_t)ldvt);
+
+                delete[] d;
+                delete[] y;
+            }
+};
+#endif

svd.h

+#ifndef SVD_H
+#define SVD_H
+
+#include "../CORE/core_types.h"
+#include "../CORE/maths.h"
+#include "../CORE/lapack.h"
+
+class SVD {
+    /** Compute Compact SVD
+     *
+     * M = USV^T
+     */
+    private:
+        char jobz = 'S'; // compute the first min(M,N) columns of U
+        int m;
+        int n;
+        int lda;
+        int ldu;
+        int ldvt;
+        int ucol;
+        int lwork;
+        int info;
+
+        // arrays
+        double *a = nullptr;
+        double *s = nullptr;
+        double *u = nullptr;
+        double *vt = nullptr;
+        double *work = nullptr;
+        int *iwork = nullptr;
+
+        bool work_alloc=false;
+
+    public:
+        ~SVD()
+        {
+            delete[] s;
+            delete[] u;
+            delete[] vt;
+
+            if (work_alloc)
+                delete[] iwork;
+
+            if (work_alloc && lwork>0)
+                delete[] work;
+        }
+        SVD(const SVD &other)
+        {
+            std::cout << "Copy-constructor" << std::endl;
+            jobz = other.jobz;
+            m = other.m;
+            n = other.n;
+            lda = other.lda;
+            ldu = other.ldu;
+            ldvt = other.ldvt;
+            ucol = other.ucol;
+            lwork = other.lwork;
+            info = other.info;
+
+            // allocate all
+            s = new double[sizeS()];
+            u = new double[sizeU()];
+            vt = new double[sizeVT()];
+
+            // We don't need this arrays
+            a = nullptr;
+            work = nullptr;
+            iwork = nullptr;
+
+            work_alloc=false;
+
+            std::cout << sizeS() << std::endl;;
+            std::cout << sizeU() << std::endl;;
+            std::cout << sizeVT() << std::endl;;
+
+            // hard copy
+            std::memcpy(s, other.s, sizeS()*sizeof(double));
+            std::memcpy(u, other.u, sizeU()*sizeof(double));
+            std::memcpy(vt, other.vt, sizeVT()*sizeof(double));
+        }
+
+        //SVD &operator=(const SVD &other)
+        //{
+        //    //if (this == &other)
+        //    //    return *this;
+        //    std::cout << "Copy-assignment" << std::endl;
+        //}
+
+        SVD(Matrix &Phi):
+            m(Phi.rows()),
+            n(Phi.cols()),
+            lda(Phi.rows()),
+            ldu(Phi.rows()),
+            ldvt(std::min(m,n)),
+            ucol(std::min(m,n)),
+            a(&Phi.data()[0]),
+            s(new double[sizeS()]),
+            u(new double[sizeU()]),
+            vt(new double[sizeVT()]),
+            iwork(new int[8*std::min(m,n)]),
+            work_alloc(true)
+        {
+            lwork=-1; // QUERRY
+            double wkopt;
+            dgesdd_(&jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
+                    &wkopt, &lwork, iwork, &info);
+            lwork = (int)wkopt;
+            work = new double[lwork];
+            dgesdd_(&jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
+                    work, &lwork, iwork, &info);
+            // Filter out v.small singularvalues from S and its inverse
+            for (size_t i=0; i<sizeS(); ++i) {
+                if (s[i]<std::numeric_limits<double>::min()) {
+                    s[i] = 0.0;
+                }
+            }
+        }
+
+        int get_info() const{
+            return info;
+        }
+        double *getA() const{
+            return a;
+        }
+        double *getU() const{
+            return u;
+        }
+        double *getVT() const{
+            return vt;
+        }
+        double *getS() const{
+            return s;
+        }
+        size_t sizeA() const {
+            return m*n;
+        }
+        size_t sizeU() const {
+            return ldu*ucol;
+        }
+        size_t sizeVT() const {
+            return ldvt*n;
+        }
+        size_t sizeS() const {
+            return std::min(m,n);
+        }
+        std::pair<size_t,size_t> shapeA() const{
+            return std::pair<size_t, size_t>(m,n);
+        }
+        std::pair<size_t,size_t> shapeU() const{
+            return std::pair<size_t, size_t>(ldu,ucol);
+        }
+        std::pair<size_t,size_t> shapeVT() const{
+            return std::pair<size_t, size_t>(ldvt,n);
+        }
+        std::pair<size_t,size_t> shapeS() const{
+            return std::pair<size_t, size_t>(std::min(m,n),1);
+        }
+
+        template <typename T>
+            T predict(const T &w) const {
+                /* Return U S V^T w */
+
+                double *t = new double[shapeA().first];
+                double *x = new double[shapeA().second];
+
+                T res = predict(w,t,x);
+
+                delete[] t;
+                delete[] x;
+                return res;
+            }
+        template <typename T>
+            T predict(const T &w, double *t, double *x) const {
+                /* Return U S V^T w */
+
+                const double *w_ptr = &w.data()[0];
+                // 1. V^T w
+                char trans='N';
+                int m = static_cast<int>(shapeA().first);
+                int n = static_cast<int>(shapeA().second);
+                double alpha_ = 1.0;
+                double beta_ = 0.0;
+                int incx=1;
+
+                dgemv_(&trans, &n, &n, &alpha_, vt, &n, w_ptr, &incx, &beta_, x, &incx);
+
+                // 2. S (V^T w)
+                for (int i=0; i<n; ++i) {
+                    x[i]*=s[i];
+                }
+
+                // 3. U (S V^T w)
+                dgemv_(&trans, &m, &n, &alpha_, u, &m, x, &incx, &beta_, t, &incx);
+                return T(t,m);
+            }
+
+};
+
+#endif