Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Linear algebra re-write #34

Merged
merged 31 commits into from
Dec 19, 2019
Merged

Linear algebra re-write #34

merged 31 commits into from
Dec 19, 2019

Conversation

mtanneau
Copy link
Member

@mtanneau mtanneau commented Dec 18, 2019

The goal of this PR is to make the linear algebra layer considerably more flexible.

Linear algebra layers

AbstractLinearSolver

  • Prior to this PR, the Newton system is systematically reduced to the normal equations, which are solved using a Cholesky factorization of the PSD matrix A*D*A'.
    In this PR, we introduce the AbstractLinearSolver type, which --among others-- makes it easier to choose between:

    • augmented system vs normal equations,
    • direct vs indirect method,
    • linear algebra backend (e.g., LAPACK, SuiteSparse, etc..).
  • Only the resolution of the augmented system is exposed to the algorithm. The reduction to the normal equations, if applicable, is done within the linear solver object itself. See the docs for more details.

  • For sparse constraint matrix (and Float32/Float64 precision), the default setting is to use an LDL factorization of the (quasi-definite) augmented system. This approach better handles dense columns in the constraint matrix, and tends to be more numerically stable.
    For dense matrices, the reduction to the normal equations is automatic, and BLAS/LAPACK are used when applicable.
    These settings are currently not exposed, i.e., one must modify the source code manually.

Regularizations

The augmented system

[ -Θ    A' ] [ x ] = [ xi_d ]
[  A    0  ] [ y ] = [ xi_p ]

is symmetric indefinite, thus its LDL factorization may not exist.
Therefore, we introduce primal and dual regularizations Rp and Rd, which make the system quasi-definite:

[ -(Θ+Rp)    A' ] [ x ] = [ xi_d ]
[      A     Rd ] [ y ] = [ xi_p ]

Here, Rp and Rd are diagonal matrices with positive diagonals, and any implementation of AbstractLinearSolver is expected to be able to handle them explicitly.

Future changes

The current interface for AbstractLinearSolver is expected to change in the near future, with the introduction of traits and extra parameters to make it more user-friendly and (even more) flexible.

Regularized algorithm

Given that we use black-box linear solvers, which does not allow explicit control over the factorization process (when a direct method is employed), the regularizations are handled in the IPM algorithm directly.

The homogeneous self-dual algorithm was modified accordingly, including updated stopping criteria. See paper for details.

* Remove commented code in HSD algorithm and related files
* Remove obsolete functions in LinearAlgebra module
* Update tests to reflect new conventions
* De-activate multiple-precision tests (temporary)
* De-activate UnitBlockAngular tests (temporary)
* Float32 is hardly used anyway
* Rational breaks type stability, since `sqrt(::Rational)` returns a `Float64`
* Computation of residuals
* Stopping criterion
* Computation of h0 in the Newton system solve
@codecov-io
Copy link

codecov-io commented Dec 18, 2019

Codecov Report

Merging #34 into master will decrease coverage by 7.75%.
The diff coverage is 81.98%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #34      +/-   ##
==========================================
- Coverage   84.91%   77.16%   -7.76%     
==========================================
  Files          20       22       +2     
  Lines        1472     1673     +201     
==========================================
+ Hits         1250     1291      +41     
- Misses        222      382     +160
Impacted Files Coverage Δ
src/LinearAlgebra/LinearAlgebra.jl 33.33% <ø> (-23.81%) ⬇️
src/Tulip.jl 100% <ø> (ø) ⬆️
src/Solvers/HSDSolver/HSDSolver.jl 90.24% <100%> (+7.7%) ⬆️
src/Solvers/Solvers.jl 75% <100%> (ø) ⬆️
src/LinearAlgebra/LinearSolvers/dense.jl 100% <100%> (ø)
src/LinearAlgebra/LinearSolvers/LinearSolvers.jl 100% <100%> (ø)
src/LinearAlgebra/LinearSolvers/sparse.jl 45.45% <45.45%> (ø)
src/Solvers/HSDSolver/hsd_step.jl 96.37% <85.71%> (-3.63%) ⬇️
src/LinearAlgebra/unitBlockAngular.jl 0% <0%> (-69.67%) ⬇️
... and 4 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update e725c42...2f62545. Read the comment docs.

* Solve proximal problem obtained from the homogeneous self-dual form
* Remove unused proximal points
* Update code for Newton system
* Remove looser criteria for infeasibility detection
* Remove log of tau and kappa
@mtanneau mtanneau merged commit 3786da7 into master Dec 19, 2019
@mtanneau mtanneau deleted the LinearSolver branch December 19, 2019 22:56
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants