diff --git a/tensorPredictors/NAMESPACE b/tensorPredictors/NAMESPACE
index 595827d..51f28dc 100644
--- a/tensorPredictors/NAMESPACE
+++ b/tensorPredictors/NAMESPACE
@@ -11,5 +11,6 @@ export(dist.subspace)
 export(matpow)
 export(matrixImage)
 export(reduce)
+export(tensor_predictor)
 import(stats)
 useDynLib(tensorPredictors, .registration = TRUE)
diff --git a/tensorPredictors/R/RMReg.R b/tensorPredictors/R/RMReg.R
index d1b40ee..d798979 100644
--- a/tensorPredictors/R/RMReg.R
+++ b/tensorPredictors/R/RMReg.R
@@ -23,7 +23,7 @@
 #' @param Z additional covariate vector (can be \code{NULL} if not required.
 #'  For regression with intercept set \code{Z = rep(1, n)})
 #' @param y univariate response vector
-#' @param lambda penalty term, if set to \code{Inf} 
+#' @param lambda penalty term, if set to \code{Inf} max lambda is computed.
 #' @param loss loss function, part of the objective function
 #' @param grad.loss gradient with respect to \eqn{B} of the loss function
 #'  (required, there is no support for numerical gradients)
@@ -64,6 +64,7 @@ RMReg <- function(X, Z, y, lambda = 0,
         Z <- matrix(0, nrow(X), 1)
         ZZiZ <- NULL
     } else {
+        if (!is.matrix(Z)) Z <- as.matrix(Z)
         # Compute (Z' Z)^{-1} Z used to solve for beta. This is constant
         # throughout and the variable name stands for "((Z' Z) Inverse) Z"
         ZZiZ <- solve(crossprod(Z, Z), t(Z))
@@ -85,7 +86,7 @@ RMReg <- function(X, Z, y, lambda = 0,
     loss1 <- loss(B1, beta, X, Z, y)
 
     # Start without, the nesterov momentum is zero anyway
-    no.nesterov <- TRUE
+    no.nesterov <- TRUE # Set to FALSE after the first iteration
     ### Repeat untill convergence
     for (iter in 1:max.iter) {
         # Extrapolation with Nesterov Momentum
diff --git a/tensorPredictors/R/matrixImage.R b/tensorPredictors/R/matrixImage.R
index 3cf4d96..ed257a2 100644
--- a/tensorPredictors/R/matrixImage.R
+++ b/tensorPredictors/R/matrixImage.R
@@ -1,28 +1,45 @@
 #' Plots a matrix as an image
 #'
 #' @param A a matrix to be plotted
+#' @param add.values boolean indicating if matrix values are to be written into
+#'  matrix element boxes
 #' @param main overall title for the plot
 #' @param sub sub-title of the plot
 #' @param interpolate a logical vector (or scalar) indicating whether to apply
 #'  linear interpolation to the image when drawing.
 #' @param ... further arguments passed to \code{\link{rasterImage}}
 #'
+#' @examples
+#' AR <- 0.5^abs(outer(1:10, 1:10, `-`))
+#' matrixImage(AR, AR > 0.2, main = "Autoregressiv Covariance")
+#'
 #' @export
-matrixImage <- function(A, main = NULL, sub = NULL, interpolate = FALSE, ...) {
-
-    # Scale values of `A` to [0, 1] with min mapped to 1 and max to 0.
-    A <- (max(A) - A) / diff(range(A))
-
+matrixImage <- function(A, add.values = FALSE,
+    main = NULL, sub = NULL, interpolate = FALSE, ...
+) {
     # plot raster image
     plot(c(0, ncol(A)), c(0, nrow(A)), type = "n", bty = "n", col = "red",
         xlab = "", ylab = "", xaxt = "n", yaxt = "n", main = main, sub = sub)
 
-    # Add X-axis giving index
-    ind <- seq(1, ncol(A), by = 1)
-    axis(1, at = ind - 0.5, labels = ind, lwd = 0, lwd.ticks = 1)
-    # Add Y-axis
-    ind <- seq(1, nrow(A))
-    axis(2, at = ind - 0.5, labels = rev(ind), lwd = 0, lwd.ticks = 1, las = 1)
+    # Scale values of `A` to [0, 1] with min mapped to 1 and max to 0.
+    S <- (max(A) - A) / diff(range(A))
 
-    rasterImage(A, 0, 0, ncol(A), nrow(A), interpolate = interpolate, ...)
+    # Raster Image ploting the matrix with element values mapped to grayscale
+    # as big elements (original matrix A) are dark and small (negative) elements
+    # are white.
+    rasterImage(S, 0, 0, ncol(A), nrow(A), interpolate = interpolate, ...)
+
+    # Add X-axis giving index
+    x <- seq(1, ncol(A), by = 1)
+    axis(1, at = x - 0.5, labels = x, lwd = 0, lwd.ticks = 1)
+    # Add Y-axis
+    y <- seq(1, nrow(A))
+    axis(2, at = y - 0.5, labels = rev(y), lwd = 0, lwd.ticks = 1, las = 1)
+
+    # Writes matrix values (in colored element grids)
+    if (any(add.values)) {
+        if (length(add.values) > 1) { A[!add.values] <- NA }
+        text(rep(x - 0.5, nrow(A)), rep(rev(y - 0.5), each = ncol(A)), A,
+             adj = 0.5, col = as.integer(S > 0.65))
+    }
 }
diff --git a/tensorPredictors/R/tensor_predictors.R b/tensorPredictors/R/tensor_predictors.R
index cad0f5e..fd0c2a7 100644
--- a/tensorPredictors/R/tensor_predictors.R
+++ b/tensorPredictors/R/tensor_predictors.R
@@ -18,6 +18,7 @@ log.likelihood <- function(par, X, Fy, Delta.inv, da, db) {
     sum(error * (error %*% Delta.inv))
 }
 
+#' @export
 tensor_predictor <- function(X, Fy, p, t, k = 1L, r = 1L, d1 = 1L, d2 = 1L,
                              method = "KPIR_LS",
                              eps1 = 1e-2, eps2 = 1e-2, maxit = 10L) {