Non-negative Ridge Regression Implementaion 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!
 A: 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.


*

*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


*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).
