I am trying to solve a singular linear least square problem:

$$minimize: \phantom{2} ||Ax-b||^2 \\ subject \phantom{2} to: \phantom{2} x \ge 0$$

Here $ A \in R^{n \times m} $, and $ n\lt m$. here m is a very large number and the matrix $A$ is dense and almost full.

The least square solution for this is equivalent to solve (maybe not equivalent due to the constraints, just to show the reasoning):

$$ A^TAx=A^Tb,\phantom{2} x \ge 0 \phantom{,,,} (*)$$

but here $A^TA$ is singular because $ n\lt m$. So my question ends up how to solve the above constrained singular linear system.

The only idea I came up with is using SVD like the following:

$$ A^TA=VSV^T=[V_1 \phantom{,,} V_2] \begin{bmatrix} S_{11} & 0 \\ 0 & 0 \\ \end{bmatrix} \begin{bmatrix} V_1^T \\ V_2^T \\ \end{bmatrix} =V_1S_{11}V_1^T$$ then setting $ x=V_1z $ leads to $ z=V_1^Tx $, and solving $(*)$ equals to solving $$S_{11}z=V_1^TA^Tb, \phantom{,} V_1z \ge0$$ which can be reformulated as, $$minimize: \phantom{2} \frac{1}{2}z^TS_{11}z-b^TAV_1z \\ subject \phantom{2} to: \phantom{2} V_1z \ge 0$$

However, unfortunately, because the matrix $A$ is large and dense, doing SVD needs large memory and is slow. When doing this using MATLAB it prompted me saying "out of memory". Is there any other more efficient way to solve this singular least square problem?

P.S. Is it possible there is no solution to the reformulated minimization problem? I know the objective function is convex since $S_{11}$ is positive definite, and of course the linear constraints are convex. But is it possible the it is infeasible, that the linear constraints might contradict with each other?

  • 2
    $\begingroup$ I don't think this should be called a singular linear system. The top-most problem is a quadratic linear-constrained program, but if $Ax=b$ actually has a solution with $x\geq0$, then it's really just a linear program. I don't think you need to do svd on $A^\top A$, have a look at en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse as described there for the under-determined, full row rank, case (the system with $A^\top A$ comes from the over-determined case). In your problem the solution isn't even necessarily unique. $\endgroup$
    – Kirill
    Mar 2 '16 at 21:05
  • $\begingroup$ Thanks for your comments. Yes, it is not a linear system. It is a least square problem with constraints. I use the linear system to explain how I came up the idea of using SVD to convert the original optimization problem to the reformulated one. . $\endgroup$
    – JW Xing
    Mar 2 '16 at 22:00
  • $\begingroup$ The reason why I want to reformulate the original optimization to the latter one is because the latter is a convex optimization but the original is not. $\endgroup$
    – JW Xing
    Mar 2 '16 at 22:02
  • 3
    $\begingroup$ I'm suspicious of your last comment: you introduce a linear change of variables, which cannot possibly turn a non-convex problem into a convex problem (linear map applied to a convex set gives a convex set, etc.). I think I see the issue: you assume that a solution for $z$ exists, i.e., the original problem attains minimum $0$ and it's feasible. If that's the case, you should just do as explained in the wikipedia article and take svd of $A$, not $A^\top A$, you'll get a linear problem without having to compute $m\times m$ svd. If not, then the second problem is not equivalent. $\endgroup$
    – Kirill
    Mar 2 '16 at 22:36
  • 4
    $\begingroup$ You won't in general be able to solve the constrained linear system of equations that you've asked about. You should start over and ask about how to solve the least squares problem with nonnegativity constraints- there are lots of approaches to that problem, but they don't use $A^{T}Ax=A^{T}b$. $\endgroup$ Mar 3 '16 at 1:22

Pretty sure that what you want is Non-negative least squares. Plenty of implementations already exist, so you shouldn't have to write any code to get going. https://en.wikipedia.org/wiki/Non-negative_least_squares


You could try my MATLAB solver PDCO which uses an interior method and will be happy that your n < m.


options.Method = 1 % the default
d1 = 1e-4
d2 = 1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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