3

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!

I-was-a-Ki
  • 33

1 Answers1

2

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

Mustafa Eisa
  • 1,322
  • 9
  • 20

  • Notice that the x^T is on the left of S since otherwise this wouldn't make sense.

    Exactly. Sorry for my mistake.

    I took the second way and it was successful! You helped me very much, and I highly appreciate your detailed and easy-to-understand explanation. Thank you so much : )

     – I-was-a-Ki Dec 19 '16 at 07:47
  • Could you provide any reference for this methods. It would be very helpful. – Manuel Dec 22 '16 at 19:24
  • You can read about projected gradient descent here: https://www.stats.ox.ac.uk/~lienart/blog_opti_pgd.html – Mustafa Eisa Dec 23 '16 at 00:45
  • Is nonnegative ridge regression not done in most straightforward way using nonnegative least squares (nnls) in which the covariate/design matrix is augmented with a diagonal matrix with the sqrt(penalties) along the diagonal and the dep variable is augmented with a nr of zeros equal to the nr of variables? As suggested in one of the answers here: https://stats.stackexchange.com/questions/203685/how-to-perform-non-negative-ridge-regression – Tom Wenseleers Jul 08 '19 at 13:52