On a Ridge regression-like problem 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.
 A: 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.
Example
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

set.seed(17)
par(mfrow=c(2,3))
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, {
    invisible(sapply(seq_len(dim(coef)[1]), 
                     function(i) 
                       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)
}
par(mfrow=c(1,1))

