Non-negative Ridge Regression Implementaion

Question

I would like x (a 100 x 1 matrix, each element must be non-negative) which minimizes the following;

$|| Mx-p ||^2 + λ xS x'$

where λ, M (a 10 x 100 matrix), p (a 10 x 1 matrix), and S (a 100 x 100 matrix) are specified.

I have already tried function "pcls" in a R package "mcgv", but I received an error message saying "Penalized model matrix must have no more columns than rows".

The implementation can be in R, C++, or anything else.
Any help/explanations would be very helpful, thank you very much!

Answer 1

First of all, I think you mean to write $$\min_x \|Mx - p\|^2_2 + \lambda x^\top S x \text{ subject to } x \geq 0$$ where the inequality constraint is evaluated component-wise. Notice that the $x^\top$ is on the left of $S$ since otherwise this wouldn't make sense.

The error message you quoted "Penalized model matrix must have no more columns than rows" is referring to the fact that $M$ has more columns (100) than rows (10); in the language of machine learning, you have more features than examples. I would take it that you simply won't be able to solve this problem with pcls the way it is. I'll provide two solutions.

1. Assuming $\lambda > 0$ and $S$ is positive definite, the problem is convex and over-determined and therefore you can solve this and guarantee a (approximate) solution using projected gradient descent. The iteration scheme looks something like this: $$x_{k+1} = \left( 2 M^\top M x_k - 2M^\top p + 2\lambda Sx_k \right)_+$$ where the function $(c)_+ := \max\{c,0\}$ is applied component-wise on vectors. This equation means that the $(k+1)^{th}$ update of the solution to the minimization problem can be obtained using the previous iterate $x_k$, some matrix multiplies, and $(\cdot)_+$ to make sure you meet your non-negativity constraint. In matlab, this would be

   x = randn(100,1) % initialize at random value
   K = 1000 % total number of iterations for projected gradient descent
   Gram = M'*M % cache the gram matrix
   Mp = M'*p % cache the inner-product matrix
   
   for i = 1:K
       x = max(2*G*x - 2*Mp + 2*lamb*S*x, 0)
   end
   

2. Another (hacky and tacky) way that might help you get around your issue and still use pcls is to simply augment your design matrix $M$ and your dependent variable vector $p$ with 90 additional rows of zeros so that the problem, at least to pcls, looks like one in which there are at least as many observations as covariates. In other words, let $$\tilde{M} = \begin{bmatrix}M \\ \mathbf{0}\end{bmatrix}\in\mathbb{R}^{100\times 100}, \ \tilde{p} = \begin{bmatrix}p \\ \mathbf{0}\end{bmatrix}\in\mathbb{R}^{100}$$ where $\mathbf{0}$ is the matrix of the appropriate size (in the first case, 90 by 100 and in the second, 90 by 1). This will lead to the same solution since $\|\tilde{M}x - \tilde{p}\|_2^2 = \|Mx - p\|_2^2$ (I'll let you work the math out on that one).