I am trying to implement a sort of Ridge regression for the following problem \begin{equation} y = a^\top X b, \end{equation} where $y\in \mathbb{R}$, $a\in\mathbb{R}^M$, $b\in\mathbb{R}^N$ and $X\in\mathbb{R}^{M\times N}$ and $a^\top$ denotes the transpose matrix of $a$.

$X$ is a matrix which coefficients are unknown, and I aim to find them by minimizing the following Lagrangian \begin{equation} \min_{X} \Vert y - a^\top X b \Vert^2 + \beta \mbox{Tr}(XX^\top) \end{equation} where $\beta > 0$ is a standard Ridge regression parameter and $\mbox{Tr}$(·) is the trace operator.

I know the regression solution without the regularization term, since it can be easily calculated by means of the Moore-Penrose inverse: \begin{equation} X = (aa^\top)^{-1}ayb^\top(bb^\top)^{-1} \end{equation}

However, I am unable to find either a closed expression for $X$ when the regularization term is included, or a proper numerical algorithm to calculate it.

Is there any reference that could shed some light on this? Thank you.

  • $\begingroup$ $aa^T$ has rank 1 and has no inverse. Similar for $b$. $\endgroup$
    – gunes
    Commented Jan 31, 2021 at 21:25
  • $\begingroup$ This is standard least squares regression, merely written in a slightly different way. So, write it out explicitly for (say) $M=N=2$ to see what's going on. If you're familiar with vector manipulations, you might recognize that $a^\prime X b = \operatorname{Tr}(ab^\prime X).$ $\endgroup$
    – whuber
    Commented Jan 31, 2021 at 21:37
  • $\begingroup$ @gunes that is true, thank you. I have to rethink that part too. $\endgroup$
    – R-bit
    Commented Jan 31, 2021 at 21:59

1 Answer 1


Writing the values for observation $i$ as $\mathbf{a}_i = (a_{i1},a_{i2},\ldots, a_{im}),$ $\mathbf{b}_i=(b_{i1},b_{i2},\ldots, b_{in})$, and $y_i$ ($i=1,2,\ldots, N$), this model can be expressed as

$$E[y_i] = \sum_{j=1}^m\sum_{k=1}^n a_{ij}X_{jk}b_{ik} = \sum_l \alpha_{il} X_l $$

where $l$ indexes the $mn$ ordered pairs $(j,k)$ and $\alpha_{il} = a_{ij}b_{ik}.$ This is a standard ordinary least squares (OLS) multiple regression problem. In statistical parlance, it consists of all the "interactions" between the $\mathbf{a}$ variables and $\mathbf{b}$ variables but has no constant or direct effects. (See the subtitle to the upper left plot in the figure below for an example of the formula.) The coefficient $X_{jk},$ to be estimated, is the interaction between component $j$ of the $\mathbf{a}$ variables and component $k$ of the $\mathbf{b}$ variables.

Moreover, although Ridge Regression is most efficiently carried out with a dedicated algorithm (to rapidly update the solutions as the relaxation parameter $\beta$ is varied), it can be executed with a standard OLS formula (as quoted in the question) after observing the Ridge Regression solution is the OLS solution with an augmented design matrix. This requires just two lines in the R code below: one to construct a formula to specify the model and the second to carry out the regression.

In tests, I find this approach to be effective primarily when $N/(mn)$ is not much greater than $1.$ Otherwise it is usually misleading, concerning both how many coefficients of $X$ to use and how to estimate them (it shrinks them too much). Consider using the Lasso instead.


I generated explanatory variables and coefficients independently with a standard Normal distribution, then added iid standard Normal noise to the computed values of the $y_i.$ I randomly zeroed out half the coefficients in $X.$

Here are ridge trace plots for two cases where $m=2,$ $n=3,$ and $N/(mn)$ is $1.5$ (top row) and $15$ (bottom row). The true coefficients (components of $X$) are plotted along the vertical axis using red for true zeros and black otherwise. (The red ticks are jittered to resolve their overlap.)


The regularization helps in the top row but seems less than useful in the bottom row.

Here is the (reproducible) R code used to implement the solution, run the simulation, and plot the results.

m <- 2
n <- 3
N <- ceiling(1.5*m*n)
sigma <- 1 # Error S.D.
p <- 1/2   # Fraction of true coefficients to be zero

for (iter in 1:6) {
  # Determine the number of observations and establish the range of 
  # coeffcients `beta` for the ridge regression.
  if (iter > 3) N <- ceiling(15*m*n)  # Increase `N` for the bottom row
  beta=seq(0, N/(m*n), length.out=101)# Specify the ridge trace (search space)
  # Generate a model `X` and data `a` and `b`.
  i <- sample.int(n*m, floor(p*n*m))            # These coefficients will be zero
  X <- rnorm(m*n); X[i] <- 0; X <- matrix(X, m) # The coefficients
  a <- matrix(rnorm(m*N, m), N)                 # The `a` explanatory variables
  b <- matrix(rnorm(n*N, n), N)                 # The `b` explanatory variables
  # Create a data frame that holds all the `a` and `b` variables.
  df <- apply(cbind(a, b), 2, scale)
  colnames(df) <- c(anames <- paste0("a", seq_len(m)), 
                    bnames <- paste0("b", seq_len(n)))
  # Compute the "true" response `y.0` and add noise for the observed response.
  y.0 <-  sapply(seq_len(N), function(i) df[i, 1:m] %*% X %*% df[i, 1:n+m])
  df <- as.data.frame(cbind(df, y=y.0 + rnorm(N, 0, sigma)))
  # Ridge regression.
  library(MASS) # lm.ridge
  # Specify the model.
  f <- as.formula(paste0("y ~ (", paste(anames, collapse="+"), ")*(",
                         paste(bnames, collapse="+"), ") - (",
                         paste(c("1", anames, bnames), collapse="+"), ")"))
  # Fit it using Ridge Regression.
  fit <- lm.ridge(f, df, lambda=beta)
  # (That's it!) ####
  # Plot the ridge trace.
  # plot(fit, xlab=expression(beta), ylab="Coefficient", lwd=2) # Built-in ridge trace
  subtitle <- gsub(" ", "", ifelse(iter==1, paste0("y ~ ", as.character(f)[3], ""), ""))
  subtitle <- sub("-", " - ", subtitle)
  plot(range(fit$lambda), range(c(c(X), fit$coef)), bty="n", type="n",
       xlab=expression(beta), ylab="Coefficient",
       main=paste0("Ridge Trace for N=", N, " and mn=", m*n), cex.main=1.2)
  mtext(subtitle, line=0.25, cex=min(0.6, 5/(m+n+1)))
  with(fit, {
                       lines(lambda, coef[i,], lwd=2, lty=(1 + i%%6), 
                                     col=hsv(i/dim(coef)[1], .8, .7))))
  # Add the rug plot.
  k <- sum(c(X)==0)
  rug(seq(-k, k, by=2)*0.01, col="Red", side=2)
  rug(c(X)[X != 0], side=2)
  • 1
    $\begingroup$ Thank you for your answer, it is really enlightening. To be honest, I had never considered the possibility of introducing the LASSO term in the objective function, but I will definitely give it a try. However, I am a bit cautious regarding its performance, due to its tendency to force some of the regression coefficients to zero. Thanks again :) $\endgroup$
    – R-bit
    Commented Feb 6, 2021 at 11:55
  • 1
    $\begingroup$ The nice thing about framing this as a standard linear regression problem is that you now have available the whole panoply of techniques, including Ridge Regression, Lasso, Elastic Net, Bayesian methods, and so on. $\endgroup$
    – whuber
    Commented Feb 6, 2021 at 17:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.