Friday, November 18, 2016

Constrained Least Squares solution times II

In this post I discussed a somewhat large linearly constrained least squares problem:

\[\begin{align}\min_P\>\>&||P s - f||^2\\ & \sum_i p_{i,j}\le 1\\ & p_{i,j} \ge 0 \end{align}\]

Note that \(i\) has 39 elements and \(j\) has 196 elements. We can solve this as a Quadratic Programming (QP) problem:

\[\begin{align}\min_x\>\>&x^TQx +c^Tx\\&Ax \{\le,=,\ge\} b\\&\ell\le x \le u\end{align}\]

There are two obvious ways to implement this. It is not immediately clear to me which one is best.

Formulation I
\[\boxed{\begin{align}\min_{P}\>&\sum_i \left(\sum_j p_{i,j} s_j – f_i \right)^2\\ & \sum_i p_{i,j} \le 1 \\ &p_{i,j} \ge 0 \end{align}}\]

For our data set this leads to a QP problem with 7644 variables and 196 linear constraints. The number of nonzero element in the lower triangle of the (quadratic) Q matrix is 401544. 

Formulation II 

\[\boxed{\begin{align}\min_{P,e}\>&\sum_i e_i^2\\ & \sum_j p_{i,j} s_j = f_i + e_i\\ &\sum_i p_{i,j} \le 1 \\ &p_{i,j} \ge 0 \end{align}}\]

Here we made the model substantially larger: we added variables \(e_i\) (the residuals) and we added (linear) constraints. Note that we created a pretty dense block here. However the Q matrix becomes much simpler: it is now a diagonal matrix. The statistics are: 7683 variables, 235 linear constraints and only 39 diagonal elements in Q.

We moved complexity from the nonlinear objective into the linear constraints. That sounds like a good thing. But we increased the number of variables, equations and nonzero-elements. The trade-off is not obvious.

Some results
Model Solver Time (sec) Notes
I Cplex 0.16 Warning:  diagonal perturbation of 1.5e-008 required to create PSD Q.
Warning:  Barrier optimality criterion affected by large objective shift.
II Cplex 0.08  
I IPOPT with MUMPS linear solver >1000 Could not solve in 1000 seconds. Messages:
MUMPS returned INFO(1) = -9 and requires more memory, reallocating.  Attempt 1
  Increasing icntl[13] from 1000 to 2000.
MUMPS returned INFO(1) = -9 and requires more memory, reallocating.  Attempt 2
  Increasing icntl[13] from 2000 to 4000.

Looks like IPOPT has some severe problems recovering after a reallocation (may be a bug).
II IPOPT with MUMPS linear solver 0.6  
I IPOPT with MA27 linear solver 200  
II IPOPT with MA27 linear solver 0.2  

The timings are seconds spend inside the solver doing real work (excluding model generation and solver reading and setup times). Indeed it looks like formulation II gives better performance and reliability. It is noted that the set \(i\) is small compared to \(j\) which makes this dataset somewhat unusual.