Start adding LP solver

[btracey]

No description provided.

[btracey]

@vladimir-ch This solution routine works almost all of the time. If pZero = 0.0 in simplex_test.go, then it passes all of the randomly generated tests. As it's currently set, to 0.9, sometimes the Phase 1 problem is returned unbounded, which should be impossible.

Sorry for the ugly code, there's a lot that needs to be fixed up.

[btracey]

If may be when an entire row of A is all zero. Investigating

[btracey]

Nope, not the problem.

[btracey]

Ah ha! Found it. Floating point error innacuracies again.

[btracey]

Well, it solves that issue, but now on to the next

[btracey]

Okay, now the problem is that some of the swaps make the basic subset of A singular (for example, that constraint isn't active in any of the sets). I started simplexSolve, which helps, but isn't yet enough and is probably the wrong thing to do anyway. I'm done for the night though.

[btracey]

@vladimir-ch I went through and cleaned up the code. The algorithm seems mostly correct. If the random problem has no zeros (see pZero in simplex_test), then the algorithm always passes. If the problem has zeros, it always passes unless it enters the bland code. Sometimes there's a linear solve failure. I'm not sure if there's a problem in the code, or if those problems really are infeasible. If they are infeasible, I'm not sure how to detect that. Any ideas?

[vladimir-ch]

As always, give me some time to look at it.

[btracey]

Of course.

[btracey]

I fixed some of the problems, and partially added some in the meantime. As usual, floating point noise is why we can't have nice things. I added a test where the checking for linear dependence fails. We need to check the condition number of the augmented matrix instead. This requires a mat64.Cond function.

[dane-unltd]

Hey, I saw that you are working on an LP solver.
I updated the predictor-corrector interior-point method I had implemented to work with the current version of the matrix package:
https://github.com/dane-unltd/optimize/tree/addlp/convex/lp
The implementation seems to be working. I also added a straightforward reference MATLAB implementation.
If you want I can push the commit to your branch.

[btracey]

Awesome! Excited to see you have an interior point solver.

I would prefer you not push it to the branch quite yet. I'm not exactly treating the commit with respect yet, as I'm still trying to debug. Once I'm pretty sure that the Simplex implementation correct, I'd love to work with you to figure out the correct API for the package.

Does your code still work if there are zeros in the A, b and C matrices? It seems like that's where the tricky issues arise.

[dane-unltd]

I think for some cases you need a modified cholesky decomposition in the interior point solver due to the (by design) ill-conditioned matrix.
I will try to see if the amount of infeasible/unbounded problems is the same as for other solvers.

[dane-unltd]

I fixed an issue in my solver regarding the step-size.
Generating problems randomly with a lot of zeros yields a lot of unbounded and infeasible problems.
I think sometimes problems are found to be infeasible which are actually unbounded.
Also when A does not have full-rank the standard cholesky decomposition does not work.

The simplex also has issues with rank-deficient A so some preprocessing is probably needed in this case.

[btracey]

Agreed on the difficulty with zeros. Simplex tries to correct for rank-deficient A, but the current code is not a sufficient detector. This is the motivation for adding Cond to mat64, but I haven't had a chance to update Simplex (been busy trying to get a lapack-based SVD into mat64).

I'm not sure if you're talking about interior point or simplex with unbounded, but you should take a look at the simplex tests. You can check the result by running the dual.

[btracey]

Okay, I think this simplex is good enough now. It passes 200,000+ random LPs with no zeros, and 50,000+ LPs when pZero = 0.7. On the 64,000 somethingth test, the condition number of ab goes bad. If the tolerance for d is decreased form 1e-13 to 1e-12, instead of the condition number going bad, instead the code never completes. Seems like the typical floating point issues.

There are a number of efficiency improvements that can happen, and I think we should find a good API for the package. (Then we can work on a QP solver, and then an SQP solver!)

[dane-unltd]

I will try to get the interior point method to pass your tests.
Regarding the API, I do not have a strong opinion.
We should avoid something like the SDPT3 matlab API, though.

I am just now implementing the modified Cholesky decomposition for the interior point method
(good overview in http://pages.cs.wisc.edu/~swright/papers/P600.pdf which also explains the interior point implementation) and will see how it goes.

[btracey]

Getting it to pass the tests isn't trivial. You can see all of the floating point number tweaks that had to go into Simplex. Hopefully interior point is more robust.

My thought was:
a) Have lp.Solve which is the easy way to call the package
b) Have lp.SolveStandardForm for those who have already paid the conversion costs
c) Export both Simplex and InteriorPoint for those that care.

I'd also like to try and modify Solve to use matrix interfaces. This way a single code could work for Dense (as it does now), but also for different Sparse implementations. Do you think IP can be similarly modified?

[dane-unltd]

Yes, I think support for sparse matrices is essential for an LP solver.
We probably want to implement fast paths for dense and sparse if a general Matrix is supplied.
Otherwise copy the provided matrix into the more convenient representation.

Of course we first need a package for sparse matrices which implements the necessary factorizations ...

[dane-unltd]

when I run the tests I get a panic for the simplex solver after i=15200
lp: unexpected linear solve error: matrix singular or near-singular with inverse condition number 6.9833e+17

I now implemented a self-dual approach from the description of the MOSEK interior point solver.
The advantage is that it can reliably detect infeasibility, at least when the matrix A is not too badly conditioned...
From what I read, they use the interior-point solution to find a basis for the problem and if it is not the optimal basis, they use it as a warm-start for the simplex method.

Extending the interior-point solvers to the sparse case would make a lot of sense but also require some work. So I am not sure if I am able to make progress in this direction any time soon.

Anyway, once you have merged this pull request, I can add another one for the interior point solvers.
Maybe we just use the API you mentioned above. It seems reasonable and could enable others to try the solvers.

[sbinet]

is this still relevant, now that we have #165 ?

[sbinet]

@btracey, can't we just close/abandon this PR?

[btracey]

@dane-unltd if you are still interested in joining forces for an interior point solver, let us know.