3
$\begingroup$

[ question reposted from https://math.stackexchange.com/questions/786612/solving-a-linearly-constrained-sparse-linear-least-squares-problem ]

Given the system of equations

$Ax=b$, subject to $Cx\le d$

where $A$ is an $n\times m$ matrix (with $n>m$) and is very large and sparse. As an example $A$ can have $3126250\times 2740$ elements. Each row of $A$ has only 4 or 5 non-zero numbers which can only be 1 or -1.

I am on Matlab and I've been using LSQR but I need the inequality constraints to impose monotonicity on $x$.

Can you please advise on any solvers to do this with linear constraints? Is there any implementation on Matlab or C for this?

$\endgroup$
  • $\begingroup$ Is your problem actually $\min \| Ax - b \|_{2}^{2}$ subject to $Cx \leq d$? $\endgroup$ – Brian Borchers May 9 '14 at 13:43
  • $\begingroup$ How big are $n$ and $m$? $\endgroup$ – Brian Borchers May 9 '14 at 13:44
  • $\begingroup$ As an example $A$ can have $3126250\times 2740$ elements. I've added more information above. My system is $Ax=b$, right now I am solving it like that yes: $argmin_x ||Ax-b||^2_2$, s.t. $Cx\le d$. $\endgroup$ – strangelyput May 9 '14 at 14:28
  • $\begingroup$ Are you constraints just of the form $x_{k} \geq \sum_{i=1}^{k-1} x_{i}$, or are there other linear inequalities? $\endgroup$ – Brian Borchers May 9 '14 at 15:22
  • $\begingroup$ They are only monotonicity constraints, $x_{i+1}\ge x_i$. $\endgroup$ – strangelyput May 12 '14 at 11:00
3
$\begingroup$

If you have access to the MATLAB optimization toolbox then this can easily be done using the quadprog() function. You'd start by writing the objective in quadratic form as

$ \| Ax - b \|_{2}^{2} = x^{T}(A^{T}A)x-2(A^{T}b)^{T}x+b^{T}b$

then multiply to get $P=A^{T}A$ and $q=-2A^{T}b$. Then your objective is

$f(x)=x^{T}Px+q^{T}x+b^{T}b$

and ready to feed into quadprog(). The $P$ matrix is only 2740 by 2740, so this isn't a very large problem from the point of view of quadprog().

I'm sure there are some free qp solvers for MATLAB if you don't have a copy of the optimization toolbox.

Also note that you may want to reformulate the problem in terms of variables $z_{i}$, where

$x_{1}=z_{1}$

$x_{2}=z_{1}+z_{2}$

$\ldots$

Then you can replace your inequality constraints $Cx \leq d$ with $z \geq 0$.

$\endgroup$
  • $\begingroup$ Thank you very much. MATLAB's quadprog with the active-set algorithm (which is the appropriate for non-convex problems with linear inequality constraints) does not support sparse matrices. I will try that approach where $z_{i+1}=x_{i+1}+x_i$ and report back. $\endgroup$ – strangelyput May 9 '14 at 18:17
  • 1
    $\begingroup$ The $A^{T}A$ matrix is small enough that it doesn't matter much whether you treat it as a sparse matrix (in fact, the matrix will likely be quite dense.) Just take full(A'*A) to get the $P$ matrix in dense form. Furthermore, the objective function will be convex, since $A^{T}A$ is positive semidefinite. $\endgroup$ – Brian Borchers May 9 '14 at 18:51
  • $\begingroup$ MATLAB's interior-point-convex for quadprog solver is returning "The problem is non-convex". So not sure what is happening there. Furthermore if that's the case, what would be the reason for quadprog to have non-convex algorithms? And yes of course you are right, there is no reason for me to use sparse with $A^T A$. Just as a note, even with a 2740 by 2740, at least with MATLAB algorithms for non-convex problems, this is prohibitively slow. $\endgroup$ – strangelyput May 12 '14 at 11:01
  • $\begingroup$ That's quite odd. I'd compute the eigenvalues of $P$ matrix next. If any of them are significantly negative then you've made a mistake in computing $P$. Its more likely that you'll find a few eigenvalues that are just barely negative (something like -1.0e-19) If that's the case, you can just add a small multiple of $I$ to $P$ to restore the positive definiteness. $\endgroup$ – Brian Borchers May 12 '14 at 13:34
  • $\begingroup$ The smaller eigenvalue is about -1e-10. I did that identity trick and it worked. Thank you so much for your help. $\endgroup$ – strangelyput May 12 '14 at 14:28

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.