# Utilizing LP constraints structure with Tulip?

**URL:** <https://discourse.julialang.org/t/utilizing-lp-constraints-structure-with-tulip/88363>\
**Category:** Optimization (Mathematical)\
**Tags:** question, linearalgebra\
**Created:** [October 6, 2022, 5:30pm UTC](https://discourse.julialang.org/t/utilizing-lp-constraints-structure-with-tulip/88363 "2022-10-06T17:30:27Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [October 6, 2022, 5:30pm UTC](https://discourse.julialang.org/t/utilizing-lp-constraints-structure-with-tulip/88363/1 "2022-10-06T17:30:27Z")

</div>

I’m [solving](https://gitlab.com/nsajko/PolynomialPassingThroughIntervals.jl) many linear optimization problems that differ in variable bounds, but share the same constraints. There is no objective, i.e, I’m just looking for a feasible solution in each case. I’m using `MultiFloats.Float64x8` arithmetic (basically a tuple of eight Float64 values).

The Tulip.jl solver seems to be able to be customized for exploiting structure in the constraint matrix, but I’d appreciate some pointers as to how to accomplish that.

This is the LP formulation notation used in Tulip’s documentation:  
 ![formulation](https://global.discourse-cdn.com/julialang/original/3X/f/4/f431a5b1a445cf68914ec817f31361b9b8dd50e7.webp)

In my case:

1. `l` and `u` (the variable bound vectors) are both finite
2. `b` is zero
3. Elements of `A` are integers, with magnitudes varying from `0` and `1` to “so big that BigInts are necessary to exactly represent it”
4. The number of variables that I’m mostly interested in is 65536 (`n == 2^16`)
5. There’s a parameter `d`, such that `m == n - d`. `d` can range from one to about thirty.
6. The matrix `A` is a block matrix formed by concatenating matrices `A0` and `I(m)`, so `A == hcat(A0, I(m))`. `I(m)` is an identity matrix of order `m`.
7. `A0` is a rectangular dense matrix with `m` rows and just `d` columns
8. Elements of `A0` are nonzero.

Are there any additional properties of my problems that I should look for?

Any advice for how to utilize this structure of my problems, if it looks like there’s some structure that might be exploitable?

---

<div class="post-metadata">

**Author:** ![ccoffrin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ccoffrin/32/400_2.png) [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)\
**Post date:** [October 7, 2022, 2:53pm UTC](https://discourse.julialang.org/t/utilizing-lp-constraints-structure-with-tulip/88363/2 "2022-10-07T14:53:02Z")

</div>

CC @mtanneau

---

<div class="post-metadata">

**Author:** ![mtanneau](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mtanneau/32/17787_2.png) [@mtanneau](https://discourse.julialang.org/u/mtanneau)\
**Post date:** [October 8, 2022, 4:03pm UTC](https://discourse.julialang.org/t/utilizing-lp-constraints-structure-with-tulip/88363/3 "2022-10-08T16:03:54Z")

</div>

Hi! Thanks for reaching out 🙂

Do you have an example of code we could run to test it out?

I see two components in your post: you’re using extended-precision arithmetic, and you’re wondering whether you can speedup the linear algebra.

Note about the first part: I’ve encountered issues with MultiFloats because of how `NaN` and `Inf` values are treated. It may not be an issue in your case because none of your bounds are infinite, but I would disable presolve just to be sure.

About the linear algebra.  
If you want to implement your custom algebra, you need the following ingredients:

- A custom type for your matrix. It will be assembled at [this line](https://github.com/ds4dm/Tulip.jl/blob/ae81ed857f96f1dbe279664f8b5504621880c61e/src/IPM/ipmdata.jl#L166):

- A custom [KKT solver](https://ds4dm.github.io/Tulip.jl/stable/reference/KKT/kkt_solvers/#AbstractKKTSolver) that can solve augmented systems, which implements the [`KKT.update!`](https://ds4dm.github.io/Tulip.jl/stable/reference/KKT/kkt_solvers/#Tulip.KKT.update!) and [`KKT.solve!`](https://ds4dm.github.io/Tulip.jl/stable/reference/KKT/kkt_solvers/#Tulip.KKT.solve!) functions.

Making a custom matrix type is the easy part. What is more involved is coding the linear solver. Namely, you need code that can solve linear systems of the form

\begin{bmatrix} -H & A^{T}\\ A & Q \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} \xi\_{d} \\ \xi\_{p} \end{bmatrix}

where H, Q are diagonal matrices with positive elements.

If I expand this system in your case, you would get:

\left[\begin{array}{ccc} -H\_{0} & & A\_{0}^{T}\\ & -H\_{s} & I \\ A\_{0} & I & Q \end{array}\right] \left[\begin{array}{c} x\_{0}\\ s\\ y \end{array} \right] = \left[\begin{array}{c} \xi^{0}\_{d}\\ \xi^{s}\_{d}\\ \xi\_{p} \end{array} \right]

where I split the variables x = (x\_{0}, s) according to the structure your described.  
The sparse linear algebra used by default will compute an LDL^{T} factorization of the above matrix, without exploiting any structure.  
In your setting, you may want to exploit the fact that A\_{0} is dense and very tall (many rows, few columns), and end up computing a dense factorization of something that looks like A\_{0}^{T}DA\_{0} + H\_{0} for some diagonal matrix D.

I [coded something similar](https://github.com/mtanneau/UnitBlockAngular.jl) for (what I called) unit block-angular matrices of the form

A = \left[\begin{array}{ccccc} B\_{0} & B\_{1} & B\_{2} &... & B\_{R}\\ 0 & e^{T} & \\ 0 && e^{T} \\ \vdots &&& \ddots\\ 0 &&&& e^{T}\\ \end{array} \right]

at the time, I observed speedups of about 10x compared to default sparse linear algebra. This was in `Float64` though, so a big chunk of speedup came from using dense BLAS calls under the hood.