Interior-point algorithms for linear programs with many inequality by bol12697

VIEWS: 7 PAGES: 12

									Interior-point algorithms for linear programs with many inequality
                            constraints

                      Luke Winternitz (lukewinternitz@gmail.com)
                                            ´
              Joint work with (advisor) Andre Tits, P.A. Absil, D.P. O’Leary




                                     April 26, 2008




 Winternitz et al. (UMCP)            IPMs w/ many constraints                  April 26, 2008   1 / 12
Primal-dual interior-point algorithms for LP


Consider the primal and dual standard LP forms
                          min c T x                    max b T y
                        s.t. Ax = b                  s.t. AT y ≤ c
                           x ≥0             (or s.t. AT y + s = c, s ≥ 0)
where A ∈ Rm×n and we assume n ≫ m.
Primal-dual interior point (PDIP) methods apply Newton’s method to the perturbed KKT
optimality conditions

                                                      AT                      c − AT y − s
                                            0                      10 1 0                  1
                                             0                   I    ∆x
       AT y + s = c,                        @A        0          0 A @∆y A = @ b − Ax A ,
               Ax = b,                       S        0          X    ∆s       σµe − Xs

              xi si = τ, i ∈ n                                             T
                                                where τ = σµ, µ = x n s > 0, σ ∈ [0, 1],
           (x, s) ≥ 0.                           X := diag(x) > 0, S := diag(s) > 0,
                                                          e := (1, 1, · · · 1)T .




     Winternitz et al. (UMCP)         IPMs w/ many constraints                      April 26, 2008   2 / 12
Cost of an iteration




Commonly, the Newton-KKT system is reduced (by block gaussian elimination) to to
the positive definite “normal-equations”

                                 AS −1 XAT ∆y = (⋆).


Dominant cost is O(nm2 ), for forming the “normal-matrix”
                                                 n
                                                 X xi
                                AS −1 XAT =           ai aT .
                                                          i
                                                   si
                                                  i=1




     Winternitz et al. (UMCP)       IPMs w/ many constraints          April 26, 2008   3 / 12
The case of n ≫ m: many inequality constraints in the dual
We expect many constraints are redundant or somehow not very relevant. We could try
to guess a good set to “pay attention to” and ignore the rest.
                      irrelevant?



                                                                                                       b
                                    b
                 ∆y
                                                                                             ∆y

         redundant
                                        active




                                                 Ignore many constraints




                                                                             max b T y
                          max b T y
                                                                                 s.t. AT y ≤ cQ
                                                                                       Q
                     s.t. AT y ≤ c
                                                                                 (AQ := [ai1 , ai2 , · · · ], ∀ ij ∈ Q)
Some prior work in 1990’s, Dantzig and Ye [1991], Tone [1993], Den Hertog [1994], but
we think there is further room to explore.

     Winternitz et al. (UMCP)                         IPMs w/ many constraints                               April 26, 2008   4 / 12
Constraint-reduced PDIP method
Given a small set Q of critical constraints, compute the PDIP direction for the reduced
primal-dual pair
                                          T
                                    min cQ xQ
                                                                            max b T y
                                 s.t. AQ xQ = b
                                                                         s.t. AT y ≤ cQ
                                                                               Q
                                     xQ ≥ 0
                   T
                  xQ sQ
i.e. let µQ =      |Q|
                          and solve

                                0 AT
                            0                      10    1 0             1
                                    Q           IQ   ∆xQ          0
                            @ AQ 0              0 A @ ∆y A = @ b −AQ xQ A ,
                               SQ 0             XQ   ∆sQ      σµQ −XQ sQ
This defines a partial search direction (∆xQ , ∆y, ∆sQ ) which is cheaper to compute
than the full PDIP direction.
                   −1
                              X xi
               AQ SQ XQ AT :=
                           Q         ai aT
                                         i    costs O(|Q|m2 ) vs. O(nm2 )
                                  si
                                          i∈Q

Issues:
     How to update (x, y, s) using this partial direction? (next slide).
     How to select Q?

      Winternitz et al. (UMCP)                    IPMs w/ many constraints                April 26, 2008   5 / 12
A possible approach: reduced normal matrix approximation (Tits et al.
[2006])



    If we solve normal equations, e.g. in affine-scaling (σ = 0), dual-feasible case:


                                 AS −1 XAT ∆y = b
                                          ∆s = −AT ∆y               (#)
                                          ∆x = −x − S −1 X ∆s.

    One possibility is to simply make the substitution:
                                         −1
                                     AQ SQ XQ AT −→ AS −1 XAT ,
                                               Q


    in (#) above, and change nothing else. Then the search direction is fully defined.

Other methods for filling out the update of (x, y, s) are possible.




      Winternitz et al. (UMCP)           IPMs w/ many constraints          April 26, 2008   6 / 12
Choosing the constraint set Q, and a convergence result




Minimal requirements on Q:
    AQ is full row rank,
    Q contains the m most nearly active constraints,

Theorem: Given this and under a nondegeneracy assumption, we have global
convergence with a local quadratic rate of two variants of the reduced PDIP method:
    Affine-scaling (Tits et al [2006])
    (Mehrotra) predictor-corrector (Winternitz et al [2007])

However, additional constraints may be needed for good numerical performance.




     Winternitz et al. (UMCP)           IPMs w/ many constraints         April 26, 2008   7 / 12
Removing a strong assumption.
Requirement 1: “AQ is full row rank” is                One approach: detect the degenerate
troublesome.                                           situation and take a ”kernel step”

                                                                               b

                                                                                      S   tep
                       a1                                                 a1       el
                                                                                 rn
                                                                               Ke

                                 a2                                                 a2




Alternatively, regularize the reduced
problem
                                      2
               max b T y − δk y                                      a1
                                                                                                Step is now well defined


                  s.t. AT y ≤ cQ
                        Q
                                                                               a2



Theorem: Under appropriate assumptions, the convergence result again holds for a
kernel-step variant, and a regularized variant that uses δk = O(µk ).


      Winternitz et al. (UMCP)            IPMs w/ many constraints                                April 26, 2008    8 / 12
Random problems
   m = 200, n = 40000
   Generate: A, b, y0 ∼ N (0, 1), s0 ∼ U(0, 1). Normalize columns of A.
   Set c = AT y0 + s0 and x0 = e.




                      CPU time                                    Iterations
                           We can evidently choose |Q| ≈ m here.
                This implies an asymptotic (large n) savings of a factor of m.


    Winternitz et al. (UMCP)           IPMs w/ many constraints                  April 26, 2008   9 / 12
Discretized semi-infinite linear programs
Semi-infinite linear programming (SILP) is a very powerful tool for engineering design.
Applications include design of, e.g., digital filters, antenna arrays, and control systems.
Can (approximately) solve SILPs by discretizing the constraint set to obtain a finite
linear program (with n ≫ m).
Example: Minimax FIR filter design.




 min t
       ˛                        ˛
  s.t. ˛H(e jω , b) − Hd (e jω )˛ ≤ t
       ˛                        ˛

     ω ∈ [0, π]




Posed as LP of size 201 × 20000. Cost to solve:

                            Standard PDIP (MPC): 34 iterations, 61 seconds.
                            Reduced PDIP (MPC): 45 iterations, 8 seconds.

      Winternitz et al. (UMCP)              IPMs w/ many constraints          April 26, 2008   10 / 12
A possible control application: LP formulation of Markov decision
processes (MDP)
A simple version of the (discounted, infinite-horizon) MDP is to solve
                                             "∞              #
                                              X k
                           min J(x0 , µ) = E     α g(xk (µ))
                                  µ
                                                         k =1

{xk } a DTFS controlled Markov chain, with states in X = {1, 2, . . . p} and controls in
U = {1, 2, . . . q} and µ is a policy (a map from states to controls).

LP formulation:                                          Using cost approx. J = Hr (r ∈ Rk ) in (∗),
                                                         get “approximate LP” (ALP) formulation:

              max e T J                                                 max (H T c)T r
             s.t. J ≤ g + αPu J       (∗)                                s.t. (Hr ) ≤ g + αPu Hr
           for all u ∈ U.                                              for all u ∈ U.

p variables and pq constraints                           k variables and pq ≫ k constaints.

Reduced PDIP methods may be useful for this class of problems.


      Winternitz et al. (UMCP)              IPMs w/ many constraints                       April 26, 2008   11 / 12
Last slide
Summary:
    We proposed a variation of the PDIP method for solving LPs with many inequality
    constraints. Basic idea is to use small subset of constraints to construct search
    direction.
    Proved global and local quadratic convergence of several variants with only
    requirement that the m most nearly active constraints are always in Q.
    Numerical results are promising.
Future work:
    Take a closer look at applications.
    Extension to convex programming, general nonlinear-programming. Conceptually,
    there is no problem if the constraints are bent outward or inward.



                                              ∇f (
                                                  y)




     Winternitz et al. (UMCP)      IPMs w/ many constraints             April 26, 2008   12 / 12

								
To top