add: cve_simple pure C implementation.

CVE_C/R/CVE.R

@@ -55,21 +55,21 @@
 #' @import stats
 #' @importFrom stats model.frame
 #' @export
-cve <- function(formula, data, method = "simple", max.dim = 10, ...) {
+cve <- function(formula, data, method = "simple", max.dim = 10L, ...) {
     # check for type of `data` if supplied and set default
     if (missing(data)) {
         data <- environment(formula)
     } else if (!is.data.frame(data)) {
-        stop('Parameter `data` must be a `data.frame` or missing.')
+        stop("Parameter 'data' must be a 'data.frame' or missing.")
     }
 
     # extract `X`, `Y` from `formula` with `data`
     model <- stats::model.frame(formula, data)
-    X <- as.matrix(model[,-1, drop = FALSE])
+    X <- as.matrix(model[ ,-1L, drop = FALSE])
-    Y <- as.matrix(model[, 1, drop = FALSE])
+    Y <- as.double(model[ , 1L])
 
     # pass extracted data on to [cve.call()]
-    dr <- cve.call(X, Y, method = method, ...)
+    dr <- cve.call(X, Y, method = method, max.dim = max.dim, ...)
 
     # overwrite `call` property from [cve.call()]
     dr$call <- match.call()
@@ -84,35 +84,81 @@ cve <- function(formula, data, method = "simple", max.dim = 10, ...) {
 #' @param ... Method specific parameters.
 #' @rdname cve
 #' @export
-cve.call <- function(X, Y, method = "simple", nObs = nrow(X)^.5,
+cve.call <- function(X, Y, method = "simple",
-                     min.dim = 1, max.dim = 10, k, ...) {
+                     nObs = sqrt(nrow(X)), h = NULL,
+                     min.dim = 1L, max.dim = 10L, k = NULL,
+                     tau = 1.0, tol = 1e-3,
+                     epochs = 50L, attempts = 10L,
+                     logger = NULL) {
 
     # parameter checking
-    if (!is.matrix(X)) {
+    if (!(is.matrix(X) && is.numeric(X))) {
-        stop('X should be a matrices.')
+        stop("Parameter 'X' should be a numeric matrices.")
     }
-    if (is.matrix(Y)) {
+    if (!is.numeric(Y)) {
-        Y <- as.vector(Y)
+        stop("Parameter 'Y' must be numeric.")
+    }
+    if (is.matrix(Y) || !is.double(Y)) {
+        Y <- as.double(Y)
     }
     if (nrow(X) != length(Y)) {
-        stop('Rows of X and number of Y elements are not compatible.')
+        stop("Rows of 'X' and 'Y' elements are not compatible.")
     }
     if (ncol(X) < 2) {
-        stop('X is one dimensional, no need for dimension reduction.')
+        stop("'X' is one dimensional, no need for dimension reduction.")
     }
 
-    if (!missing(k)) {
+    if (missing(k) || is.null(k)) {
-        min.dim <- as.integer(k)
-        max.dim <- as.integer(k)
-    } else {
         min.dim <- as.integer(min.dim)
         max.dim <- as.integer(min(max.dim, ncol(X) - 1L))
+    } else {
+        min.dim <- as.integer(k)
+        max.dim <- as.integer(k)
     }
     if (min.dim > max.dim) {
-        stop('`min.dim` bigger `max.dim`.')
+        stop("'min.dim' bigger 'max.dim'.")
     }
     if (max.dim >= ncol(X)) {
-        stop('`max.dim` must be smaller than `ncol(X)`.')
+        stop("'max.dim' (or 'k') must be smaller than 'ncol(X)'.")
+    }
+
+    if (is.function(h)) {
+        estimate.bandwidth <- h
+        h <- NULL
+    }
+
+    if (!is.numeric(tau) || length(tau) > 1L || tau <= 0.0) {
+        stop("Initial step-width 'tau' must be positive number.")
+    } else {
+        tau <- as.double(tau)
+    }
+    if (!is.numeric(tol) || length(tol) > 1L || tol < 0.0) {
+        stop("Break condition tolerance 'tol' must be not negative number.")
+    } else {
+        tol <- as.double(tol)
+    }
+
+    if (!is.numeric(epochs) || length(epochs) > 1L) {
+        stop("Parameter 'epochs' must be positive integer.")
+    } else if (!is.integer(epochs)) {
+        epochs <- as.integer(epochs)
+    }
+    if (epochs < 1L) {
+        stop("Parameter 'epochs' must be at least 1L.")
+    }
+    if (!is.numeric(attempts) || length(attempts) > 1L) {
+        stop("Parameter 'attempts' must be positive integer.")
+    } else if (!is.integer(attempts)) {
+        attempts <- as.integer(attempts)
+    }
+    if (attempts < 1L) {
+        stop("Parameter 'attempts' must be at least 1L.")
+    }
+
+    if (is.function(logger)) {
+        loggerEnv <- environment(logger)
+    } else {
+        loggerEnv <- NULL
     }
 
     # Call specified method.
@@ -120,15 +166,40 @@ cve.call <- function(X, Y, method = "simple", nObs = nrow(X)^.5,
     call <- match.call()
     dr <- list()
     for (k in min.dim:max.dim) {
+
+        if (missing(h) || is.null(h)) {
+            h <- estimate.bandwidth(X, k, nObs)
+        } else if (is.numeric(h) && h > 0.0) {
+            h <- as.double(h)
+        } else {
+            stop("Bandwidth 'h' must be positive numeric.")
+        }
+
         if (method == 'simple') {
-            dr.k <- cve_simple(X, Y, k, nObs = nObs, ...)
+            dr.k <- .Call('cve_simple', PACKAGE = 'CVE',
-        } else if (method == 'linesearch') {
+                          X, Y, k, h,
-            dr.k <- cve_linesearch(X, Y, k, nObs = nObs, ...)
+                          tau, tol,
-        } else if (method == 'sgd') {
+                          epochs, attempts,
-            dr.k <- cve_sgd(X, Y, k, nObs = nObs, ...)
+                          logger, loggerEnv)
+        #     dr.k <- cve_simple(X, Y, k, nObs = nObs, ...)
+        # } else if (method == 'linesearch') {
+        #     dr.k <- cve_linesearch(X, Y, k, nObs = nObs, ...)
+        # } else if (method == 'rcg') {
+        #     dr.k <- cve_rcg(X, Y, k, nObs = nObs, ...)
+        # } else if (method == 'momentum') {
+        #     dr.k <- cve_momentum(X, Y, k, nObs = nObs, ...)
+        # } else if (method == 'rmsprob') {
+        #     dr.k <- cve_rmsprob(X, Y, k, nObs = nObs, ...)
+        # } else if (method == 'sgdrmsprob') {
+        #     dr.k <- cve_sgdrmsprob(X, Y, k, nObs = nObs, ...)
+        # } else if (method == 'sgd') {
+        #     dr.k <- cve_sgd(X, Y, k, nObs = nObs, ...)
         } else {
             stop('Got unknown method.')
         }
+        dr.k$B <- null(dr.k$V)
+        dr.k$loss <- mean(dr.k$L)
+        dr.k$h <- h
         dr.k$k <- k
         class(dr.k) <- "cve.k"
         dr[[k]] <- dr.k
@@ -165,43 +236,20 @@ cve.call <- function(X, Y, method = "simple", nObs = nrow(X)^.5,
 #' @method plot cve
 #' @export
 plot.cve <- function(x, ...) {
+    L <- c()
+    k <- c()
+    for (dr.k in x) {
+        if (class(dr.k) == 'cve.k') {
+            k <- c(k, paste0(dr.k$k))
+            L <- c(L, dr.k$L)
+        }
+    }
+    L <- matrix(L, ncol = length(k))
+    boxplot(L, main = "Loss ...",
+            xlab = "SDR dimension k",
+            ylab = expression(L(V, X[i])),
+            names = k)
 
-
-    # H <- x$history
-    # H_1 <- H[!is.na(H[, 1]), 1]
-
-    # defaults <- list(
-    #     main = "History",
-    #     xlab = "Iterations i",
-    #     ylab = expression(loss == L[n](V^{(i)})),
-    #     xlim = c(1, nrow(H)),
-    #     ylim = c(0, max(H)),
-    #     type = "l"
-    # )
-
-    # call.plot <- match.call()
-    # keys <- names(defaults)
-    # keys <- keys[match(keys, names(call.plot)[-1], nomatch = 0) == 0]
-
-    # for (key in keys) {
-    #     call.plot[[key]] <- defaults[[key]]
-    # }
-
-    # call.plot[[1L]] <- quote(plot)
-    # call.plot$x <- quote(1:length(H_1))
-    # call.plot$y <- quote(H_1)
-
-    # eval(call.plot)
-
-    # if (ncol(H) > 1) {
-    #     for (i in 2:ncol(H)) {
-    #         H_i <- H[H[, i] > 0, i]
-    #         lines(1:length(H_i), H_i)
-    #     }
-    # }
-    # x.ends <- apply(H, 2, function(h) { length(h[!is.na(h)]) })
-    # y.ends <- apply(H, 2, function(h) { tail(h[!is.na(h)], n=1) })
-    # points(x.ends, y.ends)
 }
 
 #' Prints a summary of a \code{cve} result.

CVE_C/src/callLogger.c

@@ -0,0 +1,45 @@
+#include "cve.h"
+
+/**
+ * Calles a R function passed to the algoritm and supplied intermidiate
+ * optimization values for logging the optimization progress.
+ * The supplied parameters to the logger functions are as follows:
+ * - attempt: Attempts counter.
+ * - epoch: Current epoch staring with 0 as initial epoch.
+ * - L: Per X_i to X_j pair loss.
+ * - V: Current estimated SDR null space basis.
+ * - tau: Step-size.
+ * - err: Error \eqn{|| V V^T - V_{tau} V_{tau}^T ||}.
+ *
+ * @param logger Pointer to a SEXP R object representing an R function.
+ * @param loggerEnv Pointer to a SEXP R object representing an R environment.
+ */
+void callLogger(SEXP logger, SEXP env,
+                const int attempt, const int epoch,
+                const double* L, const int lenL,
+                const double* V, const int nrowV, const int ncolV,
+                const double tau) {
+    /* Create R objects to be passed to R logger function. */
+    // Attempt is converted from 0-indexed to 1-indexed as R index.
+    SEXP R_attempt = PROTECT(ScalarInteger(attempt + 1));
+    // No index shift for the epoch because the 0 epoch is before the first
+    // optimization step.
+    SEXP R_epoch = PROTECT(ScalarInteger(epoch));
+    SEXP R_L = PROTECT(allocVector(REALSXP, lenL));
+    SEXP R_V = PROTECT(allocMatrix(REALSXP, nrowV, ncolV));
+    SEXP R_tau = PROTECT(ScalarReal(tau));
+
+    /* Copy data to created R objects. */
+    memcpy(REAL(R_L), L, lenL * sizeof(double));
+    memcpy(REAL(R_V), V, nrowV * ncolV * sizeof(double));
+
+    /* Create logger function call as R language expression. */
+    SEXP R_exp = PROTECT(lang6(logger, R_epoch, R_attempt,
+                               R_L, R_V, R_tau));
+
+    /* Evaluate the logger function call expression. */
+    eval(R_exp, env);
+
+    /* Unprotext created R objects. */
+    UNPROTECT(6);
+}

CVE_C/src/config.h

@@ -1,7 +0,0 @@
-#ifndef CVE_INCLUDE_GUARD_CONFIG_
-#define CVE_INCLUDE_GUARD_CONFIG_
-
-#define CVE_MEM_CHUNK_SIZE 2032
-#define CVE_MEM_CHUNK_SMALL 1016
-
-#endif /* CVE_INCLUDE_GUARD_CONFIG_ */

CVE_C/src/cve.h

@@ -0,0 +1,98 @@
+/* Include Guard */
+#ifndef CVE_INCLUDE_GUARD_H_
+#define CVE_INCLUDE_GUARD_H_
+
+#include <string.h> // `mem*` functions.
+#include <math.h>   // sqrt, exp, ...
+
+#include <R.h>
+#include <Rinternals.h>
+#include <R_ext/BLAS.h>
+#include <R_ext/Lapack.h>
+
+#define CVE_MEM_CHUNK_SIZE 2032
+#define CVE_MEM_CHUNK_SMALL 1016
+
+void cve_simple_sub(const int n, const int p, const int q,
+                    const double *X, const double *Y, const double h,
+                    const double tau_init, const double tol_init,
+                    const int epochs, const int attempts,
+                    double *V, double *L,
+                    SEXP logger, SEXP loggerEnv);
+
+void callLogger(SEXP logger, SEXP env,
+                const int attempt, const int epoch,
+                const double* L, const int lenL,
+                const double* V, const int nrowV, const int ncolV,
+                const double tau);
+
+/* CVE sub-routines */
+int getWorkLen(const int n, const int p, const int q);
+double cost(const int n,
+            const double *Y,
+            const double *vecK,
+            const double *colSums,
+            double *y1, double *L);
+void scaling(const int n,
+             const double *Y, const double *y1, const double *L,
+             const double *vecD, const double *vecK,
+             const double *colSums,
+             double *vecS);
+
+/* rStiefl */
+void rStiefl(const int p, const int q, double *V,
+             double *workMem, int workLen);
+
+/* MATRIX */
+double norm(const double *A, const int nrow, const int ncol,
+            const char *type);
+
+void matrixprod(const double *A, const int nrowA, const int ncolA,
+                const double *B, const int nrowB, const int ncolB,
+                double *C);
+
+void crossprod(const double *A, const int nrowA, const int ncolA,
+               const double *B, const int nrowB, const int ncolB,
+               double *C);
+
+void skew(const int nrow, const int ncol,
+          double alpha, const double *A, const double *B,
+          double beta,
+          double *C);
+
+void nullProj(const double *B, const int nrowB, const int ncolB,
+              double *Q);
+
+void scale(const double s, double *A, const int nelem);
+
+void cayleyTransform(const int p, const int q,
+                     const double *A, const double *B,
+                     double *X, double *workMem);
+
+/* Row and column opperations. */
+void rowSums(const double *A, const int nrow, const int ncol,
+             double *sum);
+
+void colSums(const double *A, const int nrow, const int ncol,
+             double *sum);
+
+void rowSquareSums(const double *A, const int nrow, const int ncol,
+                   double *sum);
+
+void rowSumsSymVec(const double *Avec, const int nrow,
+                   const double diag,
+                   double *sum);
+
+void rowDiffs(const double* X, const int nrow, const int ncol,
+              double *diffs);
+
+void rowDiffSquareSums(const double* X, const int nrow, const int ncol,
+                       double *sum);
+
+/* SWEEP */
+void rowSweep(const double *A, const int nrow, const int ncol,
+              const char* op,
+              const double *v, // vector of length nrow
+              double *C);
+
+#endif /* CVE_INCLUDE_GUARD_H_ */

CVE_C/src/cve_simple.c

@@ -0,0 +1,165 @@
+#include "cve.h"
+
+// TODO: clarify
+static inline double gaussKernel(const double x, const double scale) {
+    return exp(scale * x * x);
+}
+
+void cve_simple_sub(const int n, const int p, const int q,
+                    const double *X, const double *Y, const double h,
+                    const double tau_init, const double tol_init,
+                    const int epochs, const int attempts,
+                    double *V, double *L,
+                    SEXP logger, SEXP loggerEnv) {
+
+    int attempt, epoch, i, j, k, nn = (n * (n - 1)) / 2;
+    double loss, loss_last, loss_best, err, tau;
+    double tol = tol_init * sqrt((double)(2 * q));
+    double gKscale = -0.5 / h;
+
+    /* Create further intermediate or internal variables. */
+    double *Q       = (double*)R_alloc(p * p,  sizeof(double));
+    double *V_best  = (double*)R_alloc(p * q,  sizeof(double));
+    double *L_best  = (double*)R_alloc(n,      sizeof(double));
+    double *V_tau   = (double*)R_alloc(p * q,  sizeof(double));
+    double *X_diff  = (double*)R_alloc(nn * p, sizeof(double));
+    double *X_proj  = (double*)R_alloc(nn * p, sizeof(double)); // TODO: needed?
+    double *y1      = (double*)R_alloc(n     , sizeof(double)); // TODO: needed?
+    double *vecD    = (double*)R_alloc(nn,     sizeof(double));
+    double *vecK    = (double*)R_alloc(nn,     sizeof(double));
+    double *vecS    = (double*)R_alloc(nn,     sizeof(double));
+    double *colSums = (double*)R_alloc(n,      sizeof(double));
+    double *G       = (double*)R_alloc(p * q,  sizeof(double));
+    double *A       = (double*)R_alloc(p * p,  sizeof(double));
+
+    /* Determine size of working memory used by subroutines. */
+    const int workLen = getWorkLen(n, p, q);
+    double *workMem = (double*)R_alloc(workLen, sizeof(double));
+
+    /* Compute X_diff, this is static for the entire algorithm. */
+    rowDiffs(X, n, p, X_diff);
+
+    for (attempt = 0; attempt < attempts; ++attempt) {
+        /* (Re)set learning rate. */
+        tau = tau_init;
+
+        /* Sample start value from stiefl manifold. */
+        rStiefl(p, q, V, workMem, workLen);
+
+        /* Create projection matrix for initial `V`. */
+        nullProj(V, p, q, Q);
+
+        /* Compute Distance vector. */
+        matrixprod(X, n, p, Q, p, p, X_proj); // here X_proj is only `(n, p)`.
+        rowDiffSquareSums(X_proj, n, p, vecD);
+
+        /* Apply kernel to distances. */
+        for (i = 0; i < nn; ++i) {
+            vecK[i] = gaussKernel(vecD[i], gKscale);
+        }
+
+        /* Compute col(row) sums of kernal vector (sym. packed lower tri
+         * matrix.), because `K == K^T` the rowSums are equal to colSums. */
+        rowSumsSymVec(vecK, n, 1.0, colSums);
+
+        /* Compute loss given the kernel vector and its column sums.
+         * Additionally the first momentum `y1` is computed and stored in
+         * the working memory (only intermidiate result, needed for `vecS`). */
+        loss_last = cost(n, Y, vecK, colSums, y1, L);
+
+        if (logger) {
+            callLogger(logger, loggerEnv,
+                       attempt, 0,
+                       L, n, V, p, q, tau);
+        }
+
+        /* Calc the scaling vector used for final computation of grad. */
+        scaling(n, Y, y1, L, vecD, vecK, colSums, vecS);
+
+        /* Compute the eucledian gradient `G`. */
+        rowSweep(X_diff, nn, p, "*", vecS, X_proj);
+        crossprod(X_diff, nn, p, X_proj, nn, p, workMem);
+        matrixprod(workMem, p, p, V, p, q, G);
+        scale(-2. / (((double)n) * h * h), G, p * q); // in-place
+
+        /* Compute Skew-Symmetric matrix `A` used in Cayley transform.
+         + `A <- tau * (G V^T - V G^T) + 0 * A`*/
+        skew(p, q, tau, G, V, 0.0, A);
+
+        for (epoch = 0; epoch < epochs; ++epoch) {
+            /* Move V allong A */
+            cayleyTransform(p, q, A, V, V_tau, workMem);
+
+            /* Create projection matrix for `V_tau`. */
+            nullProj(V_tau, p, q, Q);
+
+            /* Compute Distance vector. */
+            matrixprod(X, n, p, Q, p, p, X_proj); // here X_proj only `(n, p)`.
+            rowDiffSquareSums(X_proj, n, p, vecD);
+
+            /* Apply kernel to distances. */
+            for (i = 0; i < nn; ++i) {
+                vecK[i] = gaussKernel(vecD[i], gKscale);
+            }
+
+            /* Compute col(row) sums of kernal vector (sym. packed lower tri
+            * matrix.), because `K == K^T` the rowSums are equal to colSums. */
+            rowSumsSymVec(vecK, n, 1.0, colSums);
+
+            /* Compute loss given the kernel vector and its column sums.
+            * Additionally the first momentum `y1` is computed and stored in
+            * the working memory (only intermidiate result, needed for `vecS`). */
+            loss = cost(n, Y, vecK, colSums, y1, L);
+
+            /* Check if step is appropriate, iff not reduce learning rate. */
+            if ((loss - loss_last) > 0.0) {
+                tau *= 0.5;
+                scale(0.5, A, p * p);
+                continue;
+            }
+
+            // Compute error, use workMem (keep first `n`, they store `y1`).
+            skew(p, q, 1.0, V, V_tau, 0.0, workMem);
+            err = norm(workMem, p, p, "F");
+
+            // Shift next step to current step and store loss to last.
+            memcpy(V, V_tau, p * q * sizeof(double));
+            loss_last = loss;
+
+            if (logger) {
+                callLogger(logger, loggerEnv,
+                           attempt, epoch + 1,
+                           L, n, V, p, q, tau);
+            }
+
+            // Check Break condition.
+            if (err < tol || epoch + 1 >= epochs) {
+                break;
+            }
+
+            /* Continue computing the gradient. */
+            /* Calc the scaling vector used for final computation of grad. */
+            scaling(n, Y, y1, L, vecD, vecK, colSums, vecS);
+
+            /* Compute the eucledian gradient `G`. */
+            rowSweep(X_diff, nn, p, "*", vecS, X_proj);
+            crossprod(X_diff, nn, p, X_proj, nn, p, workMem);
+            matrixprod(workMem, p, p, V, p, q, G);
+            scale(-2. / (((double)n) * h * h), G, p * q); // in-place
+
+            /* Compute Skew-Symmetric matrix `A` used in Cayley transform.
+            + `A <- tau * (G V^T - V G^T) + 0 * A`*/
+            skew(p, q, tau, G, V, 0.0, A);
+        }
+
+        /* Check if current attempt improved previous ones */
+        if (attempt == 0 || loss < loss_best) {
+            loss_best = loss;
+            memcpy(V_best, V, p * q * sizeof(double));
+            memcpy(L_best, L, n * sizeof(double));
+        }
+    }
+
+    memcpy(V, V_best, p * q * sizeof(double));
+    memcpy(L, L_best, n * sizeof(double));
+}

CVE_C/src/cve_subroutines.c

@@ -0,0 +1,73 @@
+#include "cve.h"
+
+int getWorkLen(const int n, const int p, const int q) {
+    int mpq; /**< Max of p and q */
+    int nn = ((n - 1) * n) / 2;
+
+    if (p > q) {
+        mpq = p;
+    } else {
+        mpq = q;
+    }
+    if (nn * p < (mpq + 1) * mpq) {
+        return 2 * (mpq + 1) * mpq;
+    } else {
+        return (nn + mpq) * mpq;
+    }
+}
+
+double cost(const int n,
+            const double *Y,
+            const double *vecK,
+            const double *colSums,
+            double *y1, double *L) {
+    int i, j, k;
+    double tmp;
+
+    for (i = 0; i < n; ++i) {
+        y1[i] = Y[i];
+        L[i] = Y[i] * Y[i];
+    }
+
+    for (k = j = 0; j < n; ++j) {
+        for (i = j + 1; i < n; ++i, ++k) {
+            y1[i] += Y[j] * vecK[k];
+            y1[j] += Y[i] * vecK[k];
+            L[i] += Y[j] * Y[j] * vecK[k];
+            L[j] += Y[i] * Y[i] * vecK[k];
+        }
+    }
+
+    for (i = 0; i < n; ++i) {
+        y1[i] /= colSums[i];
+        L[i] /= colSums[i];
+    }
+
+    tmp = 0.0;
+    for (i = 0; i < n; ++i) {
+        tmp += (L[i] -= y1[i] * y1[i]);
+    }
+    return tmp / (double)n;
+}
+
+void scaling(const int n,
+             const double *Y, const double *y1, const double *L,
+             const double *vecD, const double *vecK,
+             const double *colSums,
+             double *vecS) {
+    int i, j, k, nn = (n * (n - 1)) / 2;
+    double tmp;
+
+    for (k = j = 0; j < n; ++j) {
+        for (i = j + 1; i < n; ++i, ++k) {
+            tmp = Y[j] - y1[i];
+            vecS[k]  = (L[i] - (tmp * tmp)) / colSums[i];
+            tmp = Y[i] - y1[j];
+            vecS[k] += (L[j] - (tmp * tmp)) / colSums[j];
+        }
+    }
+
+    for (k = 0; k < nn; ++k) {
+        vecS[k] *= vecK[k] * vecD[k];
+    }
+}

CVE_C/src/export.c

@@ -1,28 +1,56 @@
-#include <Rinternals.h>
+#include "cve.h"
 
-void grad(const int n, const int p, const int q,
+// SEXP rStiefl_c(SEXP pin, SEXP qin) {
-          const double *X,
+//     int p = asInteger(pin);
-          const double *X_diff,
+//     int q = asInteger(qin);
-          const double *Y,
-          const double *V,
-          const double h,
-          double *G, double *loss);
 
-SEXP grad_c(SEXP X, SEXP X_diff, SEXP Y, SEXP V, SEXP h) {
+//     SEXP Vout = PROTECT(allocMatrix(REALSXP, p, q));
-    SEXP G = PROTECT(allocMatrix(REALSXP, nrows(V), ncols(V)));
-    SEXP loss = PROTECT(ScalarReal(0.0));
 
-    grad(nrows(X), ncols(X), ncols(V),
+//     int workLen = 2 * (p + 1) * q;
-         REAL(X), REAL(X_diff), REAL(Y), REAL(V), *REAL(h),
+//     double *workMem = (double*)R_alloc(workLen, sizeof(double));
-         REAL(G), REAL(loss));
+
+//     rStiefl(p, q, REAL(Vout), workMem, workLen);
+
+//     UNPROTECT(1);
+//     return Vout;
+// }
+
+SEXP cve_simple(SEXP X, SEXP Y, SEXP k, SEXP h,
+                SEXP tau, SEXP tol,
+                SEXP epochs, SEXP attempts,
+                SEXP logger, SEXP loggerEnv) {
+    /* Handle logger parameter, set to NULL pointer if not a function. */
+    if (!(isFunction(logger) && isEnvironment(loggerEnv))) {
+        logger = (SEXP)0;
+    }
+
+    /* Get dimensions. */
+    int n = nrows(X);
+    int p = ncols(X);
+    int q = p - asInteger(k);
+
+    /* Convert types if needed. */
+    // TODO:
+
+    /* Create output list. */
+    SEXP Vout = PROTECT(allocMatrix(REALSXP, p, q));
+    SEXP Lout = PROTECT(allocVector(REALSXP, n));
+
+    /* Call CVE simple subroutine. */
+    cve_simple_sub(n, p, q,
+                   REAL(X), REAL(Y), asReal(h),
+                   asReal(tau), asReal(tol),
+                   asInteger(epochs), asInteger(attempts),
+                   REAL(Vout), REAL(Lout),
+                   logger, loggerEnv);
 
     SEXP out = PROTECT(allocVector(VECSXP, 2));
-    SET_VECTOR_ELT(out, 0, G);
+    SET_VECTOR_ELT(out, 0, Vout);
-    SET_VECTOR_ELT(out, 1, loss);
+    SET_VECTOR_ELT(out, 1, Lout);
     SEXP names = PROTECT(allocVector(STRSXP, 2));
-    SET_STRING_ELT(names, 0, mkChar("G"));
+    SET_STRING_ELT(names, 0, mkChar("V"));
-    SET_STRING_ELT(names, 1, mkChar("loss"));
+    SET_STRING_ELT(names, 1, mkChar("L"));
-    setAttrib(out, install("names"), names);
+    setAttrib(out, R_NamesSymbol, names);
 
     UNPROTECT(4);
     return out;

CVE_C/src/grad.c

@@ -1,123 +0,0 @@
-#include <stdlib.h>
-#include <math.h>
-
-#include "sums.h"
-#include "sweep.h"
-#include "matrix.h"
-#include "indexing.h"
-
-// TODO: clarify
-static inline double gaussKernel(const double x, const double scale) {
-    return exp(scale * x * x);
-}
-
-// TODO: mutch potential for optimization!!!
-static void weightedYandLoss(const int n,
-                             const double *Y,
-                             const double *vecD,
-                             const double *vecW,
-                             const double *colSums,
-                             double *y1, double *L, double *vecS,
-                             double *loss) {
-    int i, j, k, N = n * (n - 1) / 2;
-    double l;
-
-    for (i = 0; i < n; ++i) {
-        y1[i] = Y[i] / colSums[i];
-        L[i] = Y[i] * Y[i] / colSums[i];
-    }
-
-    for (k = j = 0; j < n; ++j) {
-        for (i = j + 1; i < n; ++i, ++k) {
-            y1[i] += Y[j] * vecW[k] / colSums[i];
-            y1[j] += Y[i] * vecW[k] / colSums[j];
-            L[i] += Y[j] * Y[j] * vecW[k] / colSums[i];
-            L[j] += Y[i] * Y[i] * vecW[k] / colSums[j];
-        }
-    }
-
-    l = 0.0;
-    for (i = 0; i < n; ++i) {
-        l += (L[i] -= y1[i] * y1[i]);
-    }
-    *loss = l / (double)n;
-
-    for (k = j = 0; j < n; ++j) {
-        for (i = j + 1; i < n; ++i, ++k) {
-            l = Y[j] - y1[i];
-            vecS[k]  = (L[i] - (l * l)) / colSums[i];
-            l = Y[i] - y1[j];
-            vecS[k] += (L[j] - (l * l)) / colSums[j];
-        }
-    }
-
-    for (k = 0; k < N; ++k) {
-        vecS[k] *= vecW[k] * vecD[k];
-    }
-}
-
-void grad(const int n, const int p, const int q,
-          const double *X,
-          const double *X_diff,
-          const double *Y,
-          const double *V,
-          const double h,
-          double *G, double *loss) {
-    // Number of X_i to X_j not trivial pairs.
-    int i, N = (n * (n - 1)) / 2;
-    double scale = -0.5 / h;
-
-    if (X_diff == (void*)0) {
-        // TODO: ...
-    }
-
-    // Allocate and compute projection matrix `Q = I_p - V * V^T`
-    double *Q = (double*)malloc(p * p * sizeof(double));
-    nullProj(V, p, q, Q);
-
-    // allocate and compute vectorized distance matrix with a temporary
-    // projection of `X_diff`.
-    double *vecD = (double*)malloc(N * sizeof(double));
-    double *X_proj;
-    if (p < 5) { // TODO: refine that!
-        X_proj = (double*)malloc(N * 5 * sizeof(double));
-    } else {
-        X_proj = (double*)malloc(N * p * sizeof(double));
-    }
-    matrixprod(X_diff, N, p, Q, p, p, X_proj);
-    rowSquareSums(X_proj, N, p, vecD);
-
-    // Apply kernel to distence vector for weights computation.
-    double *vecW = X_proj; // reuse memory area, no longer needed.
-    for (i = 0; i < N; ++i) {
-        vecW[i] = gaussKernel(vecD[i], scale);
-    }
-    double *colSums = X_proj + N; // still allocated!
-    rowSumsSymVec(vecW, n, 1.0, colSums); // rowSums = colSums cause Sym
-
-    // compute weighted responces of first end second momontum, aka y1, y2.
-    double *y1 = X_proj + N + n;
-    double *L = X_proj + N + (2 * n);
-    // Allocate X_diff scaling vector `vecS`, not in `X_proj` mem area because
-    // used symultanious to X_proj in final gradient computation.
-    double *vecS = (double*)malloc(N * sizeof(double));
-    weightedYandLoss(n, Y, vecD, vecW, colSums, y1, L, vecS, loss);
-
-    // compute the gradient using X_proj for intermidiate scaled X_diff.
-    rowSweep(X_diff, N, p, "*", vecS, X_proj);
-    // reuse Q which has the required dim (p, p).
-    crossprod(X_diff, N, p, X_proj, N, p, Q);
-    // Product with V
-    matrixprod(Q, p, p, V, p, q, G);
-    // And final scaling (TODO: move into matrixprod!)
-    scale = -2.0 / (((double)n) * h * h);
-    N = p * q;
-    for (i = 0; i < N; ++i) {
-        G[i] *= scale;
-    }
-
-    free(vecS);
-    free(X_proj);
-    free(vecD);
-    free(Q);
-}

CVE_C/src/indexing.c

@@ -1,12 +0,0 @@
-#include "indexing.h"
-
-void rangePairs(const int from, const int to, int *pairs) {
-    int i, j, k;
-
-    for (k = 0, i = from; i < to; ++i) {
-       for (j = i + 1; j < to; ++j, k += 2) {
-            pairs[k] = i;
-            pairs[k + 1] = j;
-        }
-    }
-}

CVE_C/src/indexing.h

@@ -1,8 +0,0 @@
-
-/* Include Guard */
-#ifndef CVE_INCLUDE_GUARD_INDEXING_
-#define CVE_INCLUDE_GUARD_INDEXING_
-
-void rangePairs(const int from, const int to, int *pairs);
-
-#endif /* CVE_INCLUDE_GUARD_INDEXING_ */

CVE_C/src/init.c

@@ -8,10 +8,14 @@
 */
 
 /* .Call calls */
-extern SEXP grad_c(SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP cve_simple(SEXP X, SEXP Y, SEXP k,
+                       SEXP h,
+                       SEXP tau, SEXP tol,
+                       SEXP epochs, SEXP attempts,
+                       SEXP logger, SEXP loggerEnv);
 
 static const R_CallMethodDef CallEntries[] = {
-    {"grad_c",          (DL_FUNC) &grad_c,          5},
+    {"cve_simple", (DL_FUNC) &cve_simple, 10},
     {NULL, NULL, 0}
 };
 

CVE_C/src/matrix.c

@@ -1,9 +1,26 @@
-#include <string.h> // for `mem*` functions.
+#include "cve.h"
 
-#include "config.h"
+double norm(const double *A, const int nrow, const int ncol,
-#include "matrix.h"
+            const char *type) {
+    int i, nelem = nrow * ncol;
+    int nelemb = (nelem / 4) * 4;
+    double sum = 0.0;
 
-#include <R_ext/BLAS.h>
+    if (*type == 'F') {
+        for (i = 0; i < nelemb; i += 4) {
+            sum += A[i]     * A[i]
+                 + A[i + 1] * A[i + 1]
+                 + A[i + 2] * A[i + 2]
+                 + A[i + 3] * A[i + 3];
+        }
+        for (; i < nelem; ++i) {
+            sum += A[i] * A[i];
+        }
+        return sqrt(sum);
+    } else {
+        error("Unknown norm type.");
+    }
+}
 
 void matrixprod(const double *A, const int nrowA, const int ncolA,
                 const double *B, const int nrowB, const int ncolB,
@@ -35,12 +52,12 @@ void nullProj(const double *B, const int nrowB, const int ncolB,
               double *Q) {
     const double minusOne = -1.0;
     const double one = 1.0;
+    int i, nelem = nrowB * nrowB;
 
     // Initialize as identity matrix.
     memset(Q, 0, sizeof(double) * nrowB * nrowB);
-    double *Q_diag, *Q_end = Q + nrowB * nrowB;
+    for (i = 0; i < nelem; i += nrowB + 1) {
-    for (Q_diag = Q; Q_diag < Q_end; Q_diag += nrowB + 1) {
+        Q[i] = 1.0;
-        *Q_diag = 1.0;
     }
 
     // DGEMM with parameterization:
@@ -50,10 +67,29 @@ void nullProj(const double *B, const int nrowB, const int ncolB,
                     &one, Q, &nrowB);
 }
 
+// In place scaling of elements of A.
+void scale(const double s, double *A, const int nelem) {
+    int i, nelemb = (nelem / 4) * 4;
+
+    if (s == 1.) {
+        return;
+    }
+
+    for (i = 0; i < nelemb; i += 4) {
+        A[i]     *= s;
+        A[i + 1] *= s;
+        A[i + 2] *= s;
+        A[i + 3] *= s;
+    }
+    for (; i < nelem; ++i) {
+        A[i] *= s;
+    }
+}
+
 // A dence skwe-symmetric rank 2 update.
 // Perform the update
 //      C := alpha (A * B^T - B * A^T) + beta C
-void skewSymRank2k(const int nrow, const int ncol,
+void skew(const int nrow, const int ncol,
           double alpha, const double *A, const double *B,
           double beta,
           double *C) {
@@ -69,3 +105,60 @@ void skewSymRank2k(const int nrow, const int ncol,
              &alpha, B, &nrow, A, &nrow,
              &beta, C, &nrow);
 }
+
+
+/** Cayley transformation of matrix `B` using the Skew-Symmetri matrix `A`.
+ *      X = (I + A)^-1 (I - A) B
+ * by solving the following linear equation:
+ *      (I + A) X = (I - A) B ==> X = (I + A)^-1 (I - A) B
+ *      \_____/     \_____/
+ *        IpA   X =   ImA   B
+ *                  \_______/
+ *        IpA   X =     Y     ==>  X = IpA^-1 Y
+ *
+ * @param A Skew-Symmetric matrix of dimension `(n, n)`.
+ * @param B Matrix of dimensions `(n, m)` with `m <= n`.
+ * @return Transformed matrix `X`.
+ * @note This opperation is equivalent to the R expression:
+ *      solve(diag(1, n) + A) %*% (diag(1, n) - A) %*% B
+ * or
+ *      solve(diag(1, n) + A, (diag(1, n) - A) %*% B)
+ */
+void cayleyTransform(const int p, const int q,
+                     const double *A, const double *B,
+                     double *X, double *workMem) {
+    int i, info, pp = p * p;
+    double zero = 0., one = 1.;
+
+    /* Allocate row permutation array used by `dgesv` */
+    int *ipiv = (int*)workMem;
+
+    /* Create Matrix IpA = I + A (I plus A) */
+    double *IpA = (double*)(ipiv + p);
+    memcpy(IpA, A, pp * sizeof(double));
+    for (i = 0; i < pp; i += p + 1) {
+        IpA[i] += 1.; // +1 to diagonal elements.
+    }
+
+    /* Create Matrix ImA = I - A (I minus A) */
+    double *ImA = IpA + pp;
+    for (i = 0; i < pp; ++i) {
+        ImA[i] = -A[i];
+    }
+    for (i = 0; i < pp; i += p + 1) {
+        ImA[i] += 1.; // +1 to diagonal elements.
+    }
+
+    /* Y as matrix-matrix product of ImA and B:
+     * Y = 1 * ImA * B + 0 * Y */
+    F77_CALL(dgemm)("N", "N", &p, &q, &p,
+                    &one, ImA, &p, B, &p, &zero, X, &p);
+
+    /* Solve system IpA Y = X for Y (and store result in X).
+     * aka. X = IpA^-1 X */
+    F77_CALL(dgesv)(&p, &q, IpA, &p, ipiv, X, &p, &info);
+
+    if (info) {
+        error("[ cayleyTransform ] error in dgesv - info %d", info);
+    }
+}

CVE_C/src/matrix.h

@@ -1,25 +0,0 @@
-
-/* Include Guard */
-#ifndef CVE_INCLUDE_GUARD_MATRIX_
-#define CVE_INCLUDE_GUARD_MATRIX_
-
-void matrixprod(const double *A, const int nrowA, const int ncolA,
-                const double *B, const int nrowB, const int ncolB,
-                       double *C);
-
-void crossprod(const double *A, const int nrowA, const int ncolA,
-               const double *B, const int nrowB, const int ncolB,
-               double *C);
-
-void nullProj(const double *B, const int nrowB, const int ncolB,
-              double *Q);
-
-// A dence skwe-symmetric rank 2 update.
-// Perform the update
-//      C := alpha (A * B^T - B * A^T) + beta C
-void skewSymRank2k(const int nrow, const int ncol,
-                   double alpha, const double *A, const double *B,
-                   double beta,
-                   double *C);
-
-#endif /* CVE_INCLUDE_GUARD_MATRIX_ */

CVE_C/src/rStiefl.c

@@ -0,0 +1,81 @@
+#include "cve.h"
+
+// /**
+//  * Performas a QR factorization and computes the Q factor.
+//  * 
+//  * @param A matrix.
+//  * @returns The Q factor of the QR factorization `A = QR`.
+// */
+// SEXP qrQ(SEXP Ain) {
+//     int i, j, info;
+
+//     if (!isMatrix(Ain)) {
+//         error("Argument must be a matrix.");
+//     }
+//     int nrow = nrows(Ain);
+//     int ncol = ncols(Ain);
+
+//     double *A = (double*)R_alloc(nrow * ncol, sizeof(double));
+//     memcpy(A, REAL(Ain), nrow * ncol * sizeof(double));
+
+//     // double *A = REAL(Ain);
+//     // Scalar factors of elementary reflectors.
+//     double *tau = (double*)R_alloc(ncol, sizeof(double));
+
+//     // Create Working memory area.
+//     int lenWork = 3 * nrow;
+//     double *memWork = (double*)R_alloc(lenWork, sizeof(double));
+
+//     F77_NAME(dgeqrf)(&nrow, &ncol, A, &nrow, tau,
+//                      memWork, &lenWork, &info);
+
+//     SEXP Qout = PROTECT(allocMatrix(REALSXP, nrow, ncol));
+//     double *Q = REAL(Qout);
+
+//     for (j = 0; j < ncol; ++j) {
+//         for (i = 0; i < nrow; ++i) {
+//             if (i == j) {
+//                 Q[i + nrow * j] = 1.;
+//             } else {
+//                 Q[i + nrow * j] = 0.;
+//             }
+//         }
+//     }
+
+//     F77_NAME(dormqr)("L", "N", &nrow, &ncol, &ncol, A, &nrow, tau, Q, &nrow,
+//                      memWork, &lenWork, &info);
+
+//     UNPROTECT(1);
+//     return Qout;
+// }
+
+void rStiefl(const int p, const int q, double *V,
+             double *workMem, int workLen) {
+    int i, j, info;
+    int pq = p * q;
+
+    GetRNGstate();
+    for (i = 0; i < pq; ++i) {
+        workMem[i] = norm_rand();
+    }
+    PutRNGstate();
+
+    double *tau = workMem + pq;
+    workLen -= pq + q;
+
+    F77_CALL(dgeqrf)(&p, &q, workMem, &p, tau,
+                     workMem + pq + q, &workLen, &info);
+
+    for (j = 0; j < q; ++j) {
+        for (i = 0; i < p; ++i) {
+            if (i == j) {
+                V[i + p * j] = 1.;
+            } else {
+                V[i + p * j] = 0.;
+            }
+        }
+    }
+
+    F77_NAME(dormqr)("L", "N", &p, &q, &q, workMem, &p, tau, V, &p,
+                     workMem + pq + q, &workLen, &info);
+}

CVE_C/src/rowColOp.c

@@ -1,8 +1,12 @@
-#include <string.h> // for `mem*` functions.
+#include "cve.h"
-
-#include "config.h"
-#include "sums.h"
 
+/**
+ * Computes the row sums of a matrix `A`.
+ * @param A Pointer to col-major matrix elements, size is `nrow * ncol`.
+ * @param nrow Number of rows of `A`.
+ * @param ncol Number of columns of `A`.
+ * @param sum Pointer to output row sums of size `nrow`.
+ */
 void rowSums(const double *A, const int nrow, const int ncol,
              double *sum) {
     int i, j, block_size, block_size_i;
@@ -41,29 +45,38 @@ void rowSums(const double *A, const int nrow, const int ncol,
 }
 
 void colSums(const double *A, const int nrow, const int ncol,
-             double *sum) {
+             double *colSums) {
-    int j;
+    int i, j;
-    double *sum_end = sum + ncol;
+    int nrowb = 4 * (nrow / 4); // 4 * floor(nrow / 4)
+    double colSum;
 
-    memset(sum, 0, sizeof(double) * ncol);
+    for (j = 0; j < ncol; ++j) {
-    for (; sum < sum_end; ++sum) {
+        colSum = 0.0;
-        for (j = 0; j < nrow; ++j) {
+        for (i = 0; i < nrowb; i += 4) {
-            *sum += A[j];
+            colSum += A[i]
+                    + A[i + 1]
+                    + A[i + 2]
+                    + A[i + 3];
         }
+        for (; i < nrow; ++i) {
+            colSum += A[i];
+        }
+        *(colSums++) = colSum;
         A += nrow;
     }
 }
 
-void rowSquareSums(const double *A, const int nrow, const int ncol,
+void rowSquareSums(const double *A,
+                   const int nrow, const int ncol,
                    double *sum) {
     int i, j, block_size, block_size_i;
     const double *A_block = A;
     const double *A_end = A + nrow * ncol;
 
-    if (nrow < CVE_MEM_CHUNK_SIZE) {
+    if (nrow > CVE_MEM_CHUNK_SIZE) {
-        block_size = nrow;
-    } else {
         block_size = CVE_MEM_CHUNK_SIZE;
+    } else {
+        block_size = nrow;
     }
 
     // Iterate `(block_size_i, ncol)` submatrix blocks.
@@ -111,3 +124,37 @@ void rowSumsSymVec(const double *Avec, const int nrow,
         }
     }
 }
+
+void rowDiffs(const double* X, const int nrow, const int ncol,
+              double *diffs) {
+    int i, j, k, l;
+    const double *Xcol;
+
+    for (k = l = 0; l < ncol; ++l) {
+        Xcol = X + l * nrow;
+        for (i = 0; i < nrow; ++i) {
+            for (j = i + 1; j < nrow; ++j) {
+                diffs[k++] = Xcol[i] - Xcol[j];
+            }
+        }
+    }
+}
+
+void rowDiffSquareSums(const double* X, const int nrow, const int ncol,
+                       double *sum) {
+    int i, j, k, l;
+    const double *Xcol;
+    double tmp;
+
+    memset(sum, 0, ((nrow * (nrow - 1)) / 2) * sizeof(double));
+
+    for (l = 0; l < ncol; ++l) {
+        Xcol = X + l * nrow;
+        for (k = i = 0; i < nrow; ++i) {
+            for (j = i + 1; j < nrow; ++j, ++k) {
+                tmp = Xcol[i] - Xcol[j];
+                sum[k] += tmp * tmp;
+            }
+        }
+    }
+}

CVE_C/src/sums.h

@@ -1,19 +0,0 @@
-
-/* Include Guard */
-#ifndef CVE_INCLUDE_GUARD_SUMS_
-#define CVE_INCLUDE_GUARD_SUMS_
-
-void rowSums(const double *A, const int nrow, const int ncol,
-             double *sum);
-
-void colSums(const double *A, const int nrow, const int ncol,
-             double *sum);
-
-void rowSquareSums(const double *A, const int nrow, const int ncol,
-                   double *sum);
-
-void rowSumsSymVec(const double *Avec, const int nrow,
-                   const double diag,
-                   double *sum);
-
-#endif /* CVE_INCLUDE_GUARD_SUMS_ */

CVE_C/src/sweep.c

@@ -1,7 +1,4 @@
-#include <R_ext/Error.h> // for `error`.
+#include "cve.h"
-
-#include "config.h"
-#include "sweep.h"
 
 /* C[, j] = A[, j] * v for each j = 1 to ncol */
 void rowSweep(const double *A, const int nrow, const int ncol,
@@ -107,7 +104,5 @@ void rowSweep(const double *A, const int nrow, const int ncol,
             C_block += block_size_i;
             v += block_size_i;
         }
-    } else {
-        error("Got unknown 'op' (opperation) argument.");
     }
 }

CVE_C/src/sweep.h

@@ -1,11 +0,0 @@
-
-/* Include Guard */
-#ifndef CVE_INCLUDE_GUARD_SWEEP_
-#define CVE_INCLUDE_GUARD_SWEEP_
-
-void rowSweep(const double *A, const int nrow, const int ncol,
-              const char* op,
-              const double *v, // vector of length nrow
-              double *C);
-
-#endif /* CVE_INCLUDE_GUARD_SWEEP_ */

test.R

@@ -1,23 +1,35 @@
+
+args <- commandArgs(TRUE)
+if (length(args) > 0) {
+    method <- args[1]
+} else {
+    method <- "simple"
+}
+epochs <- 50L
+attempts <- 25L
+
 # library(CVEpureR)
-# path <- '~/Projects/CVE/tmp/logger.R.pdf'
+# path <- paste0('~/Projects/CVE/tmp/logger_', method, '.R.pdf')
-
 library(CVE)
-path <- '~/Projects/CVE/tmp/logger.C.pdf'
+path <- paste0('~/Projects/CVE/tmp/logger_', method, '.C.pdf')
 
-epochs <- 100
+# Define logger for `cve()` method.
-attempts <- 25
+logger <- function(epoch, attempt, L, V, tau) {
-
-# Define the logger for the `cve()` method.
-logger <- function(env) {
     # Note the `<<-` assignement!
-    loss.history[env$epoch + 1, env$attempt] <<- env$loss
+    loss.history[epoch + 1, attempt] <<- mean(L)
-    error.history[env$epoch + 1, env$attempt] <<- env$error
+    if (epoch == 0) {
-    tau.history[env$epoch + 1, env$attempt] <<- env$tau
+        error <- NA
+    } else {
+        error <- norm(V %*% t(V) - V_last %*% t(V_last), type = 'F')
+    }
+    V_last <<- V
+    error.history[epoch + 1, attempt] <<- error
+    tau.history[epoch + 1, attempt] <<- tau
     # Compute true error by comparing to the true `B`
-    B.est <- null(env$V) # Function provided by CVE
+    B.est <- null(V) # Function provided by CVE
     P.est <- B.est %*% solve(t(B.est) %*% B.est) %*% t(B.est)
     true.error <- norm(P - P.est, 'F') / sqrt(2 * k)
-    true.error.history[env$epoch + 1, env$attempt] <<- true.error
+    true.error.history[epoch + 1, attempt] <<- true.error
 }
 
 pdf(path)
@@ -37,12 +49,15 @@ for (name in paste0("M", seq(5))) {
     P <- B %*% solve(t(B) %*% B) %*% t(B)
 
     # Setup histories.
+    V_last <- NULL
     loss.history       <- matrix(NA, epochs + 1, attempts)
     error.history      <- matrix(NA, epochs + 1, attempts)
     tau.history        <- matrix(NA, epochs + 1, attempts)
     true.error.history <- matrix(NA, epochs + 1, attempts)
 
-    dr <- cve(Y ~ X, k = k, logger = logger, epochs = epochs, attempts = attempts)
+    dr <- cve(Y ~ X, k = k, method = method,
+              epochs = epochs, attempts = attempts,
+              logger = logger)
 
     # Plot history's
     matplot(loss.history,       type = 'l', log = 'y', xlab = 'i (iteration)',
@@ -58,3 +73,5 @@ for (name in paste0("M", seq(5))) {
         main = paste('learning rate', name),
         ylab = expression(tau[i]))
 }
+
+cat("Created plot:", path, "\n")