Linear programming with matrix constraints

Question

I have an optimization problem that looks like the following

$$ \begin{array}{rl} \min_{J,B} & \sum_{ij} |J_{ij}|\\ \textrm{s.t.} & MJ + BY =X \end{array} $$

Here, my variables are matrices $J$ and $B$, but the entire problem is still a linear program; the remaining variables are fixed.

When I attempt to enter this program into my favorite linear programming tools, I run into some trouble. Namely, if I write this in "standard" linear program form, the parameter matrices $M$ and $Y$ end up getting repeated a ton of times (once for each column of $X$).

Is there an algorithm and/or package that can deal with optimizations of the form above? Right now I'm running out of memory because $M$ and $Y$ have to be copied so many times!

Answer 1

Overview

You might want to try a variant of the Alternating Directions Method of Multipliers (ADMM), which has been found to converge surprisingly quickly for $l_1$ lasso type problems. The strategy is to formulate the problem with an augmented Lagrangian and then do gradient ascent on the dual problem. It is especially nice for this particular $l^1$ regularized problem because the nonsmooth part of each iteration of the method has an exact solution you can simply evaluate element by element, while the smooth part involves solving a linear system.

In this post we

• derive an overall ADMM formulation for a generalization of your problem,
• derive the subproblems for each ADMM iteration and specialize them to your situation, and then
• investigate the resulting linear system that needs to be solved each iteration, and develop a fast solver (or preconditioner) based on precomputing the eigenvalue decompositions (or low rank approximations thereof) for $M^TM$ and $YY^T$.
• summarize with a few concluding remarks

Most of the big ideas here are covered in the following superb review paper,

Boyd, Stephen, et al. "Distributed optimization and statistical learning via the alternating direction method of multipliers." Foundations and Trends® in Machine Learning 3.1 (2011): 1-122. http://www.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf

Before going into details, I want to note that this is a method/algorithm answer not a practical existing code answer - if you want to use this method, you would need to roll your own implementation.

ADMM formulation

In general, suppose you want to solve $$ \begin{array}{rl} \min_{x} & \sum_{i} |x_i|\\ \textrm{s.t.} & Ax = b \end{array}. $$

The problem in the original post falls in this category after appropriate vectorization. (this is only in principle - we will see that the vectorization doesn't need to be performed in practice)

You could instead solve the equivalent problem, $$ \begin{array}{rl} \min_{x,z} & \sum_{i} |x_i| + \frac{\alpha}{2}||x-z||^2 + \frac{\beta}{2}||Az-b||^2 \\ \textrm{s.t.} & Az = b \\ \textrm{&} & x = z, \end{array} $$ which has Lagrangian $$\begin{align} L(x,z,\lambda,\gamma) =& \sum_{i} |x_i| + \frac{\alpha}{2}||x-z||^2 + \frac{\beta}{2}||Az-b||^2 + \lambda^T(Az-b) + \gamma^T(x-z) \\ =& \sum_{i} |x_i| + \frac{\alpha}{2}||x-z + \frac{1}{\alpha}\gamma||^2 + \frac{\beta}{2}||Az-b + \frac{1}{\beta}\lambda||^2 \\ &+ \frac{\alpha}{2}||\frac{1}{\alpha}\gamma||^2 + \frac{\beta}{2}||\frac{1}{\beta}\lambda||^2. \end{align}$$

The alternating direction method of multipliers solves the dual problem, $$\max_{\lambda,\gamma} \min_{x,z} L(x,z,\lambda,\gamma),$$ via gradient ascent on the dual variables, except with inexact alternating projections on the dual subproblems. Ie., one does the iteration $$ \begin{align} x^{k+1} &= \mathrm{argmin}_x L(x,z^k,\lambda^k,\gamma^k) \\ z^{k+1} &= \mathrm{argmin}_z L(x^{k+1},z,\lambda^k,\gamma^k) \\ \gamma^{k+1} &= \gamma^k + \alpha(x^{k+1}-z^{k+1}) \\ \lambda^{k+1} &= \lambda^k + \beta(Az^{k+1}-b). \end{align}$$

Under certain mild conditions on the parameters $\alpha$ and $\beta$ (explained in the Boyd & Parikh paper linked above), the ADMM method will converge to the true solution. The convergence rate is linear, as it is at the core a gradient ascent method. Often it can be accelerated to be superlinear by 1) changing the parameters $\alpha$ and $\beta$ as you go based on heuristics, or 2) using Nesterov acceleration. For notes on changing the penalty parameters, see the Boyd survey paper, and for using Nesterov acceleration with ADMM see the following paper,

Goldstein, Tom, Brendan O’Donoghue, and Simon Setzer. "Fast alternating direction optimization methods." CAM report (2012): 12-35. ftp://ftp.math.ucla.edu/pub/camreport/cam12-35.pdf

However, even if the overall convergence rate is only linear, for $l^1$ problems the method has been observed to find the sparsity pattern very quickly, and then converge more slowly on the exact values. Since finding the sparsity pattern is the hardest part, this is very fortuitous! The exact reasons why seem to be an area of current research. Everyone sees the sparsity pattern converge fast, but no one seems to know exactly why it happens. A while ago I asked Boyd and Parikh about this over email and Parikh thought it might be explained by interpreting the method in a control systems context. Another heuristic explanation of the phenomenon is found in the appendix of the following paper,

Goldstein, Tom, and Stanley Osher. "The split Bregman method for L1-regularized problems." SIAM Journal on Imaging Sciences 2.2 (2009): 323-343. ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf

Of course now the difficulty lies in solving the $x$ and $z$ update subproblems for your specific situation. Since the Lagrangian is quadratic in $z$, the $z$ update subproblem simply requires solving a linear system. The $x$ subproblem seems harder since it is nondifferentiable, but it turns out that there is an exact formula for the solution that can be applied element by element! We now discuss these subproblems in more detail and specify them to the problem in the original post.

Setup for the $z$ update subproblem (linear system)

For the $z$ update, we have $$\mathrm{argmin}_z L(x_k,z,\lambda_k,\gamma_k) = \mathrm{argmin}_z \frac{\alpha}{2}||x-z + \frac{1}{\alpha}\gamma||^2 + \frac{\beta}{2}||Az-b + \frac{1}{\beta}\lambda||^2.$$

Specialized to your problem this becomes, $$\begin{align} \mathrm{argmin}_{Z_J,Z_B} &\frac{\alpha}{2}||J^{k+1}-Z_J + \frac{1}{\alpha}\Gamma_J||_{Fro}^2 + \frac{\alpha}{2}||B^{k+1}-Z_B + \frac{1}{\alpha}\Gamma_B||_{Fro}^2 \\ &+\frac{\beta}{2}||MZ_J + Z_BY - X + \frac{1}{\alpha}\Lambda||^2_{Fro}, \end{align}$$

where $||\cdot||Fro$ denotes the Frobenius (elementwise $l_2$) norm. This is a quadratic minimization problem, where the first order optimality conditions can be found by taking partial derivatives of the objective with respect to $Z_J$ and $Z_B$ and setting them to zero. This is, $$\begin{align} 0 &= -\frac{\alpha}{2}(J^{k+1} - Z_J + \frac{1}{\alpha}\Gamma_J) + \frac{\beta}{2}M^T(MZ_J + Z_BY - X + \frac{1}{\beta}\Lambda), \\ 0 &= -\frac{\alpha}{2}(B^{k+1} - Z_B + \frac{1}{\alpha}\Gamma_B) + \frac{\beta}{2}(MZ_J + Z_BY - X + \frac{1}{\beta}\Lambda)Y^T. \end{align}$$

As noted in the comments by the original poster Justin Solomon, this system for $Z_J,Z_B$ is symmetric so conjugate gradient is an ideal matrix-free method. A later section discusses this system and how to solve/precondition it in greater detail.

Solving $x$ update subproblem (analytic thresholding solution)

Now we turn to the $x$ subproblem, $$\mathrm{argmin}_x L(x,z^k,\lambda^k,\gamma^k) = \mathrm{argmin}_x \sum_{i} |x_i| + \frac{\alpha}{2}||x-z^k + \frac{1}{\alpha}\gamma^k||^2$$

The first thing to see is that the sum can be broken up element by element, $$\sum_{i} |x_i| + \frac{\alpha}{2}||x-z^k + \frac{1}{\alpha}\gamma^k||^2 = \sum_{i} |x_i| + \frac{\alpha}{2}\sum_i (x_i-z_i^k + \frac{1}{\alpha}\gamma_i^k)^2,$$

So we can solve the optimization problem element by element in parallel, yielding $$x_i^{k+1} = \mathrm{argmin}_{x_i} |x_i| + \frac{\alpha}{2}(x_i-z_i^k + \frac{1}{\alpha}\gamma_i^k)^2.$$

The general form of this equation is, $$\min_s |s| + \frac{\alpha}{2}(s-t)^2.$$

The absolute value function is trying to pull the optimal point towards $s=0$, whereas the quadratic term is trying to pull the optimal point towards $s=t$. the true solution therefore lies somewhere on the segment $[0,t)$ between the two, with increasing $\alpha$ tending to pull the optimal point towards $t$, and decreasing $\alpha$ pulling the optimal point towards $0$.

This is a convex function but it is not differentiable at zero. The condition for a minimizing point is that the subderivative of the objective at that point contains zero. The quadratic term has derivative $\alpha(s-t)$, and the absolute value function has derivative $-1$ for $s < 0$, set-valued subderivative as the interval $[-1,1]$ when $s=0$, and derivative $1$ for $s > 0$. Thus we get the subderivative for the overall objective function, $$\partial_s \left(|s| + \frac{\alpha}{2}(s-t)^2\right) = \begin{cases} 1 + \alpha (s-t)\, & s > 0 \\ [-1,1] + \alpha t, & s = 0, \\ -1 + \alpha (s-t), & s < 0. \end{cases}$$

From this we see that the subderivative of the objective at $s=0$ contains $0$ if and only if $|t| \le \frac{1}{\alpha}$, in which case $s=0$ is the minimizer. On the other hand, if $s=0$ is not the minimizer, then we can set the single-valued derivative equal to zero and solve for the minimizer. Doing this yields, $$\mathrm{argmin}_s |s| + \frac{\alpha}{2}(s-t)^2 = \begin{cases} t - \frac{1}{\alpha}, & t > \frac{1}{\alpha}, \\ 0, & |t| \le \frac{1}{\alpha}, \\ t + \frac{1}{\alpha}, & t < -\frac{1}{\alpha} \end{cases}$$

Specializing this result again to the real problem we are trying to solve in the original question where $t = Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k$ yields, $$J_{ij}^{k+1} = \begin{cases} Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k - \frac{1}{\alpha}, & Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k > \frac{1}{\alpha}, \\ 0, & |Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k| \le \frac{1}{\alpha}, \\ Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k + \frac{1}{\alpha}, & Z_{ij}^k - \frac{1}{\alpha}\Gamma_{ij}^k < -\frac{1}{\alpha}. \end{cases}$$ The update for $B$ is simply $$B^{k+1} = Z_B - \frac{1}{\alpha}\Gamma_B, $$

as noted by original poster Justin Solomon in the comments. Overall, doing the update for $J,B$ just requires looping through the entries of your matrices and evaluating the above formulas for each entry.

Schur complement for the $Z_J,Z_B$ system

The most costly step of the iteration is solving the system, $$\begin{align} 0 &= -\frac{\alpha}{2}(J^{k+1} - Z_J + \frac{1}{\alpha}\Gamma_J) + \frac{\beta}{2}M^T(MZ_J + Z_BY - X + \frac{1}{\beta}\Lambda), \\ 0 &= -\frac{\alpha}{2}(B^{k+1} - Z_B + \frac{1}{\alpha}\Gamma_B) + \frac{\beta}{2}(MZ_J + Z_BY - X + \frac{1}{\beta}\Lambda)Y^T. \end{align}$$

To that end, it is worth some effort to construct a good solver/preconditioner for this system. In this section we do so by vectorizing, forming a Schur complement, doing some Krnoecker product manipulations and then unvectorizing. The resulting Schur complement system is a slightly modified Sylvester equation.

In what follows the following identities about vectorization and Kronecker products are absolutely key:

• $\mathrm{vec}(ABC) = (C^T \otimes A)\mathrm{vec}(B),$
• $(A \otimes B)(C \otimes D) = AC \otimes BD$,
• $(A \otimes B)^{-1} = A^{-1} \otimes B^{-1}$, and
• $(A \otimes B)^T = A^T \otimes B^T$.

These identities hold whenever the matrix sizes and invertibility are such that each side of the equation is a valid expression.

The vectorized form of the system is, $$\left(\alpha I +\beta\begin{bmatrix}I \otimes M^TM & (Y \otimes M)^T \\ Y \otimes M & YY^T \otimes I\end{bmatrix}\right)\begin{bmatrix}\mathrm{vec}(Z_J) \\ \mathrm{vec}(Z_B)\end{bmatrix} = \begin{bmatrix}\mathrm{vec}(\alpha J + \beta M^TX + \Gamma_J - M^T\Lambda) \\ \mathrm{vec}(\alpha B + \beta XY^T + \Gamma_B - \Lambda Y^T)\end{bmatrix},$$

or, $$\begin{bmatrix}I \otimes (\alpha I + \beta M^TM) & \beta (Y \otimes M)^T \\ \beta Y \otimes M & (\alpha I + \beta YY^T) \otimes I\end{bmatrix} \begin{bmatrix}\mathrm{vec}(Z_J) \\ \mathrm{vec}(Z_B)\end{bmatrix} = \begin{bmatrix}\mathrm{vec}(F) \\ \mathrm{vec}(G)\end{bmatrix},$$

where $F$ and $G$ are condensed notation for the right-hand-side. Now we perform block-gaussian-elimination/Schur complement to eliminate the lower left block of the matrix, in the process condensing the Kronecker products. This is, $$\begin{bmatrix}I \otimes (\alpha I + \beta M^TM) & \beta (Y \otimes M)^T \\ 0 & (\alpha I + \beta YY^T) \otimes I - \beta^2 YY^T \otimes M(\alpha I + \beta M^TM)^{-1} M^T\end{bmatrix} \dots \\ \cdot \begin{bmatrix}\mathrm{vec}(Z_J) \\ \mathrm{vec}(Z_B)\end{bmatrix} = \begin{bmatrix}\mathrm{vec}(F) \\ \mathrm{vec}(G) - \beta Y \otimes M(\alpha I + \beta M^TM)^{-1}\mathrm{vec}(F)\end{bmatrix}.$$

Unvectorizing, the two equations we have to solve in sequence are,

1. $$Z_B (\alpha I + \beta YY^T) - (\beta M (\alpha I + \beta M^TM)^{-1} M^T)Z_B(\beta YY^T) \dots \\ = G - \beta M (\alpha I + \beta M^TM)^{-1} F Y^T$$
2. $$(\alpha I + \beta M^TM) Z_J = F - \beta M^T Z_B Y.$$

Solution of Schur complement system when $Y,M$ are squareish, high rank

In this section we solve the Schur complement system for $Z_B$ (equation 1. above) by using precomputed full SVD's of the matrices $YY^T, MM^T, M^TM$ and applying a modified version of the Bartels-Stewart algorithm for the Sylvester equation. The algorithm is slightly modified from the standard version to account for the extra $\beta YY^T$ on the second term, which makes it not quite the Sylvester equation. Once $Z_B$ is found via the first equation, $Z_J$ can be found from the second equation easily. The second equation is trivial to solve via any method you like.

This method requires an upfront cost to precompute two full SVDs before the ADMM process starts, but then is fast to apply in the actual ADMM iterations. Since the method deals with full SVDs of the constraint matrices, it is appropriate when they are close to square and high rank. A more complicated method using low rank SVD's is also possible, but is presented in a later section.

The method develops as follows. Let $$Q D Q^T = YY^T, \\ W\Sigma W^T = MM^T, \\ VTV^T = M^TM$$ denote precomputed full singular value decompositions, and condense the right hand side to be $H$. Then the first equation becomes, $$Z_B Q (\alpha I + D) Q^T - W \beta \Sigma (\alpha I + \Sigma)^{-1}\Sigma W^T Z_B Q D Q^T = H.$$ Multiplying by the orthogonal factors to clear out the left and right and setting a new temporary unknown $A = W^T Z_B Q$, this further becomes, $$A (\alpha I + D) - \beta \Sigma (\alpha I + \Sigma)^{-1}\Sigma A D = W H Q^T.$$

Now we can find $A$ by solving the diagonal system, $$\left((\alpha I + D) \otimes I + D \otimes \beta \Sigma (\alpha I + \Sigma)^{-1}\Sigma \right)\mathrm{vec}(A) = \mathrm{vec}(W H Q^T).$$

Having found $A$, we compute $Z_B = W A Q^T$, and knowing $Z_B$ we solve the second equation above for $Z_J$, which is trivial since we already have the eigenvalue decomposition for $M^TM$.

The upfront cost is computing two symmetric positive definite eigenvalue decompositions of $M^TM$ and $YY^T$, and then the cost-per-iteration for a complete solve is dominated by a handful of matrix-matrix multiplications, which is on the same order of magnitude as doing 1 CG subiteration. If the upfront eigenvalue decompositions are too costly, then they can be computed inexactly by, for example, terminating the Lanczos iteration early and keeping the largest eigenvectors. Then the method can be used as a good preconditioner for CG rather than a direct solver.

Solution method when $M,Y$ are very rectangular or have low rank approximation

Now we turn our attention to solving or preconditioning the $Z_J,Z_B$ when either a) the input matrices $M,Y$ are very rectangular - meaning they have many more rows than columns or vice versa - or b) they have low rank approximation. The derivation below involves extensive use of Woodbury formula, Schur complement, and other similar manipulations.

We start with our Schur complement system, $$(\alpha I + \beta YY^T) \otimes I - \beta^2 YY^T \otimes M(\alpha I + \beta M^TM)^{-1} M^T.$$

A few manipulations transform this system into a more symmetric form, $$(\alpha I + \beta I \otimes MM^T + \beta YY^T \otimes I)\mathrm{vec}(Z_B) = \left(I \otimes (I + \frac{\beta}{\alpha}MM^T)\right)\mathrm{vec}(H).$$

Now we bring in the low rank approximations. Let $$Q D^{1/2} Q_2^T = Y \\ W \Sigma^{1/2} V^T = M$$ be either the reduced SVD's or low rank approximations of $Y$ and $M$ ($Q_2$ is a placeholder and is not used). Substituting these into our system yields the following matrix inverse we wish to apply, $$(\alpha I + \beta I \otimes W \Sigma W^T + \beta YY^T \otimes I)^{-1}.$$

Since the matrix we with to invert is a low-rank update to the identity, the logical strategy is to try to use the Woodbury formula, $$(A + UCU^T)^{-1} = A^{-1} - A^{-1}U(C^{-1}+U^TA^{-1}U)^{-1}U^TA^{-1}.$$

However, some care is needed since the low rank pieces $I \otimes W$ and $Y \otimes I$ are not orthogonal. Thus to apply the Woodbury formula we collect both low rank updates into a single big update. Doint so and applying the Woodbury formula yields, $$\left(\frac{1}{\alpha} I + \beta \begin{bmatrix}I\otimes W & Q \otimes I\end{bmatrix}\begin{bmatrix}I \otimes \Sigma & \\ & D \otimes Y\end{bmatrix}\begin{bmatrix}I \otimes \Sigma^T \\ Q^T \otimes I\end{bmatrix}\right)^{-1} \\ = \alpha I - \frac{\beta}{\alpha^2}\begin{bmatrix}I\otimes W & Q \otimes I\end{bmatrix}\begin{bmatrix}I \otimes (\Sigma^{-1}+\frac{\beta}{\alpha}I) & \frac{\beta}{\alpha}Q \otimes W^T\\ \frac{\beta}{\alpha}Q^T\otimes W & (D^{-1} + \frac{\beta}{\alpha}I) \otimes Y\end{bmatrix}^{-1}\begin{bmatrix}I \otimes \Sigma^T \\ Q^T \otimes I\end{bmatrix}.$$

The core inverse can be computed by the blockwise 2x2 inverse formula, $$\begin{bmatrix}A & B \\ B^T & C\end{bmatrix}^{-1} = \begin{bmatrix}(A-BC^{-1}B^T)^{-1} & -A^{-1}B(C-B^TA^{-1}B)^{-1} \\ -C^{-1}B^T(A-BC^{-1}B^T)^{-1} & (C-B^TA^{-1}B)^{-1}\end{bmatrix}.$$

This post is already long enough so I'll spare the lengthy details of the calculation, but the end result is that plugging the necessary submatrices into the blockwise inverse and multiplying everything through yields the following explicit form for the overall inverse, $$(\alpha I + \beta I \otimes MM^T + \beta YY^T \otimes I)^{-1} = \frac{1}{\alpha} I - \frac{\beta}{\alpha^2}(t_{11} + s_{11} + t_{12} + s_{12} + t_{21} + s_{21} + t_{22} + s_{22}),$$

where $$\begin{align} t_{11} &= \frac{\alpha}{\beta}I \otimes W l^{-1} W^T \\ s_{11} &= (Q \otimes W l^{-1})D_{11}(Q^T \otimes l^{-1}W^T) \\ t_{12} &= -\frac{\alpha}{\beta} Q h^{-1} Q^T \otimes W l^{-1} W^T \\ s_{12} &= -(Q h^{-1} \otimes W l^{-1})D_{22}(h^{-1} Q^T \otimes W^T) \\ t_{21} &= t_{12} \\ s_{21} &= -(Q h^{-1} \otimes W)D_{22}(h^{-1} Q^T \otimes l^{-1} W^T) \\ t_{22} &= \frac{\alpha}{\beta}Q h^{-1} Q^T \otimes I \\ s_{22} &= (Q h^{-1} \otimes W)D_{22}(h^{-1}Q^T \otimes W^T) \\ D_{11} &= \frac{\alpha}{\beta}\left(h \otimes I - I \otimes l^{-1} \right)^{-1} \\ D_{22} &= \frac{\alpha}{\beta}\left(I \otimes l - h^{-1} \otimes I \right)^{-1} \\ l &= \frac{\alpha}{\beta} \Sigma^{-1} + I \\ h &= \frac{\alpha}{\beta} D^{-1} + I. \end{align}$$

In this form, we can apply the inverse and find $Z_B$ term by term through 8 left and right matrix multiplication sandwiches. The general formula for applying the sum of Kronecker products is, $$\left((A_1 \otimes B_1) + (A_2 \otimes B_2) + \dots\right)\mathrm{vec}(C) = \mathrm{vec}(B_1^T C A_1 + B_2^T C A_2 + \dots ).$$

Note that all the explicit inverses we ended out with are diagonal, so there is nothing to be "solved".

Linear solver code

I implemented the above two $z_J,Z_B$ solvers in Matlab. The seem to work well. The solver code is here.

https://github.com/NickAlger/MeshADMM/blob/master/zkronsolve.m

A test script for checking that the solvers work is here. It also shows by example how calll the solver code.

https://github.com/NickAlger/MeshADMM/blob/master/test_zkronsolve.m

Concluding remarks

ADMM-type methods are well suited for problems like this, but you would need to roll your own implementation. The overall structure of the method is pretty simple so implementation is not too difficult in something like MATLAB.

The piece missing from this post that would be need to be specified to fully define the method for your problem is a choice of penalty parameters $\alpha,\beta$. Luckily the method is generally pretty robust as long as the parameter vales aren't crazy. The Boyd and Parikh paper has a section about the penalty parameters as do the references therein, but I would just experiment with the parameters until you get reasonable convergence rates.

The $Z_J,Z_B$ solver strategies presented are highly effective if the constraint matrices are either a) dense, squareish, and high rank, or b) have a good low rank approximation. Another useful solver that could be a topic of future work would be a solver optimized for the following case - the constraint matrix $M$ is sparse and squareish and high rank, but there exists a good preconditioner for $\alpha I + MM^T$. This would be the case if, for example, $M$ is a discretized Laplacian.

Answer 2

You probably want to use a matrix-free method for linear programming. I don't know of any method specifically geared towards linear programming, but there exist matrix-free interior point methods for quadratic programs and for general nonlinear programs. The quadratic program case corresponds exactly to your problem, where the quadratic form coefficients are all zeroes. (You could also tailor methods that use exact linear solves to the structure of your problem, but that sort of bespoke implementation may not be worth it, and is less practical than using a matrix-free method.)

I don't know of any commercial optimization package that implements matrix-free variants of interior point methods. IPOPT is supposed to implement a matrix-free interior point method for nonlinear programming, but I haven't been able to track down the API calls that enable you to use it.

In addition to CVX, you could probably use GAMS or AMPL to input the matrix once and set up your constraints in the modeling language to use that matrix. However, the methods used by the solver backends to CVX, GAMS, and AMPL do not use matrix-free solvers; all will require the full coefficient matrix for the linear program in standard form, which will be huge (it will be a Kronecker product of matrices). What will probably happen is that you input your linear program in the form above using the modeling language, and then the modeling language will translate the data into a form usable by backend solvers. This form will require huge matrices, and I suspect you will run into the same sorts of errors (unless you run on a machine with enough memory).

Answer 3

Can you afford those SVDs Geoffrey Irving mentioned? If you can, I would consider an iteratively reweighted least squares (IRLS) approach. This approach would solve problems of the form $$\begin{array}{ll}\text{minimize}&\sum_{ij} W_{ij}J_{ij}^2\\\text{subject to}&MJ+BY=X\end{array}$$where $W$ is a weight matrix.

The iterations begin with $W^{(0)}$ as the all ones matrix; this yields an optimum $J^{(0)}$. The iterations proceed with $$W_{ij}^{(k+1)}=\left|\max\{J_{ij}^{(k)},\epsilon\}\right|^{-1}$$ where $\epsilon$ is a small constant that prevents division by zero. I'm not totally sure about the convergence criteria, but perhaps the Wikipedia link I offered above can give you references.

You could also consider a smoothed first-order method. TFOCS, which I co-authored, could handle this using its "smoothed conic dual" (SCD) solver, but it's not going to be as easy to use.

If you do want to try a matrix-free interior point method, read the work of Jacek Gondzio.

EDIT: hmm, it might be the case that IRLS won't be able to use the SVD to compute solutions. If so I would fall back to one of the other choices.

Answer 4

You could try using CVX, which would allow you to code it in exactly the form you've written it (i.e., with $X$ as a matrix rather than a vector). It probably would be solved using a more general convex solver instead of an LP solver, but if the convex solver succeeds then I suppose you won't mind.

Another possibility: use a solver that allows your constraint matrices to be stored as sparse matrices. This will still require a lot more memory than what you ought to need, but far less than if you stored them as dense matrices. In CVX, if you use kron you get a sparse matrix, so it would be trivial to try this.