Histórico matemático para PDLP

Esta página contém conhecimentos matemáticos para desenvolvedores e usuários avançados do PDLP, um solucionador de programação linear e quadrático disponível em ferramentas OR. Ele serve como referência para partes do código e não precisa ser lido por conta própria. Os leitores interessados precisam primeiro se familiarizar com o artigo "Programação linear prática em grande escala usando gradiente híbrido Primal-Dual e analisar o código e voltar a este documento quando o código o referenciar.

Primal

O PDLP considera o seguinte problema de programação quadrática convexa:

min

onde A é uma matriz m \times n e Q é uma matriz n \times n não negativa e uma diagonal1. Os vetores do limite superior u^{c} e u^{v} têm entradas em \mathbb{R} \cup \{ \infty\} e os vetores do limite inferior l^{c} e l^{v} têm entradas em \mathbb{R} \cup \{le} -\infty\^ ^$.

Duplo

Permita que a \in \mathbb{R}. Permita que [a]_+ indique a parte positiva e [a]_- anote a parte negativa, ou seja, a = [a]_+ - [a]_-. Quando aplicadas a um vetor, as partes positivas e negativas são calculadas com base em elementos.

O duplo do problema primário anterior está acima de x \in \mathbb{R}^n, y \in \mathbb{R}^m e r \in \mathbb{R}^n. O vetor y contém multiplicadores duplos nas restrições lineares (l^{c} \le Ax \in \mathbb{R}^m, r \in \mathbb{R}^n. O vetor y contém multiplicadores duplos nas restrições lineares (l^{c} \le Ax \in \mathbb{R}^m) contém \u$

\begin{align} \max_{x, y, r} & \, -\frac{1}{2}x^TQx + \left((l^{c})^T[y]_+ - (u^{c})^T[y]_- \right) + \left((l^{v})^T[r]_+ - (u^{v})^T[r]_- \right) \\ \text{s.t.}\; & Qx + c - A^Ty = r \end{align}

Quando Q = 0, x pode ser descartado da dupla, recuperando a dualidade de LP.

Limites de variável dupla

Dizemos que y satisfaz os limites de variáveis duplas se o y-term no objetivo for finito, ou seja:

y_i \geq 0 \qquad \text{if }u^c_i = \infty, \\ y_i \leq 0 \qquad \text{if }l^c_i = -\infty.

Derivação usando a dualidade conjugada

Informações preliminares

Deixe a \in \mathbb{R} \cup \{-\infty\} e b \in \mathbb{R} \cup \{\infty\} com b \ge a e considere o intervalo $[a, b] \subseteq \mathbb{R} \cup \in $.infty,

Permita que \mathcal{I}_{[a, b]} : \mathbb{R} \to \mathbb{R} \cup \{ \infty\} seja a função indicador do intervalo, ou seja, \mathcal{I}_{[a, b]}(x) é zero quando x \in [a, b]

Defina p(y; a, b): \mathbb{R} \to \mathbb{R} \cup \{\infty\} como:

p(y; a, b) = \begin{cases} ay & \text{ if } y < 0 \\ 0 & \text{ if } y = 0 \\ by & \text{ if } y > 0 \end{cases}.

Quando a ou b são infinitos, siga a aritmética real estendida padrão.

Resultado básico: p(y; a, b) = (\mathcal{I}_{[a, b]})^*(y), em que (\cdot)^* indica o conjugado convexo.

For vectors l \subseteq (\mathbb{R} \cup \{-\infty\})^n and u \subseteq (\mathbb{R} \cup \{\infty\})^n, the indicator function \mathcal{I}_{[l, u]} : \mathbb{R}^n \to \mathbb{R} \cup \{ \infty\} is defined analogously as \mathcal{I}_{[l,u]}(x) = \sum_{i=1}^n \mathcal{I}_{[l_i, u_i]}(x_i), and similarly p(y; l, u) = (\mathcal{I}_{[l, u]})^*(y) = \sum_{i=1}^n p(y_i; l_i, u_i).

Derivação

Introduzindo as variáveis auxiliares \tilde a \in \mathbb{R}^m e \tilde x \in \mathbb{R}^n, reafirmamos o problema primal como:

\begin{align} \min_{x, \tilde x, \tilde a} & \, c^Tx + \frac{1}{2}x^TQx + \mathcal{I}_{[l^c,u^c]}(\tilde a) + \mathcal{I}_{[l^v,u^v]}(\tilde x) \\ \text{s.t.}\; & \tilde a = Ax \\ & \tilde x = x \end{align}

Ao duplicar as restrições de igualdade, temos:

\min_{x, \tilde x, \tilde a} \max_{y, r} c^Tx + \frac{1}{2}x^TQx + y^T\tilde a - y^TAx + r^T\tilde x - r^Tx + \mathcal{I}_{[l^c,u^c]}(\tilde a) + \mathcal{I}_{[l^v,u^v]}(\tilde x)

Troca de mínimo por máximo e reagrupamento:

\max_{y, r} \min_{x, \tilde x, \tilde a} c^Tx + \frac{1}{2}x^TQx - y^TAx - r^Tx + \left( \mathcal{I}_{[l^c,u^c]}(\tilde a) + y^T\tilde a \right) + \left(\mathcal{I}_{[l^v,u^v]}(\tilde x) + r^T\tilde x \right)

A minimização de articulações acima de x, \tilde x e \tilde a se decompõe. Para x, vemos que um minimizador, se existir, satisfaz Qx + c - A^Ty = r. Nesse caso, o valor mínimo é -\frac{1}{2} x^TQx. Para \tilde x e \tilde a, aplicamos a definição de conjugados convexos com pequenos ajustes para sinais.

Isso produz o seguinte:

\begin{align} \max_{x, y, r} & \, -\frac{1}{2}x^TQx - p(-y, l^c, u^c) - p(-r, l^v, u^v) \\ \text{s.t.}\; & Qx + c - A^Ty = r \end{align}

Ao expandir a definição de p, obtemos o duplo indicado no topo.

Formulação de sela

O gradiente híbrido primitiva (consulte Chambolle e Pock) visa um problema primário do formato

\begin{align} \min_x f(x) + g(x) + h(Kx) \end{align}

que, pela dualidade conjugada, é equivalente ao problema do ponto de sela

\begin{align} \min_x \max_y f(x) + g(x) + y^TKx - h^*(y) \end{align}

O PDLP força o problema de programação quadrática convexa a este formulário ao definir:

  • f(x) = 0
  • g(x) = c^T x + \frac{1}{2} x^T Q x + \mathcal{I}_{[l^v, u^v]}(x)
  • h(a) = \mathcal{I}_{[-u^c,-l^c]}(a)
  • K = -A

Como derivado anteriormente, h^*(y) = p(y; -u^c,-l^c) é uma função convex linear em partes. Tanto g quanto h^* podem aceitar valores infinitos, o que limita efetivamente os domínios de x e y, respectivamente.

Observe que o operador proximal de g pode ser calculado em formato fechado, segundo a suposição do PDLP de que Q é diagonal, usando o fato de que g é separável e a seguinte propriedade válida para qualquer função f_1, f_2:

f_2(t) = f_1(t) + \frac{\mu}{2} t^2 \Rightarrow \mathrm{prox}_{\lambda f_2}(t) = \mathrm{prox}_{\frac{\lambda}{1 + \lambda \mu} f_1}\left( \frac{t}{1 + \lambda \mu} \right).

Para comprovar esse fato, veja o exemplo do Teorema 6.13 em Métodos de primeira ordem na otimização (link em inglês). A expressão resultante é dada

\begin{equation} \mathrm{prox}_{\lambda g}(x) = \mathrm{proj}_{[l^v, u^v]}\left( (I + \lambda Q)^{-1} (x - \lambda c) \right) \end{equation}

Custos reduzidos, resíduos duplos e a correção do objetivo duplo

A formulação de pontos de sela funciona explicitamente apenas com (x,y). Os custos reduzidos r são implícitos. Para retornar (x,y,r) ao resolver a formulação de ponto sela, definimos r como r = Qx + c - A^Ty. O objetivo duplo corrigido é o valor objetivo do problema duplo e sempre fornece um limite inferior no valor do objetivo, mas é -\infty sempre que houver um custo não infinito reduzido.

Os custos reduzidos e os objetivos duplos informados pelo PDLP foram alterados entre as versões 9.7 e 9.8 das ferramentas OR. Para verificar qual versão se aplica, verifique o comentário que descreve SolverResult::reduced_costs em primal_dual_hybrid_gradient.h e veja se ela menciona a versão 9.8.

Versão 9.7 e anteriores

Para ter um valor duplo mais significativo quando o objetivo duplo corrigido é -\infty, também informamos um objetivo duplo que ignora os termos infinitos no valor do objetivo. Os residuais duplos são os valores de r de termos infinitos no objetivo duplo corrigido, com 0 nos outros componentes, e os custos reduzidos retornados pelo PDLP são os valores de r dos termos finitos no objetivo duplo corrigido, com 0 nos outros componentes, de modo que $r = \mbox{residuals} + \mbox{d$}.

Versão 9.8 e mais recentes

Para ter um valor duplo mais significativo quando o objetivo duplo corrigido for -\infty, também informamos um objetivo duplo que substitui os termos infinitos no valor do objetivo por finitos da seguinte maneira: se um dos limites for finito, esse limite será usado no lugar do infinito. Caso contrário, o zero será usado para o limite. Essa escolha preserva a concavidade do objetivo duplo, mas não fornece necessariamente um limite inferior para o valor do objetivo. Os resíduos duais são os valores de r de termos infinitos no objetivo duplo corrigido. Os custos reduzidos retornados pelo PDLP são de r.

Tratamento de alguns limites variáveis como infinito

Nas duas versões, se a opção handle_some_primal_gradients_on_finite_bounds_as_residuals do solucionador for true (o padrão), limites de variáveis adicionais poderão ser tratados como infinitos ao calcular o objetivo duplo e os resíduos duplos. Em particular, se |x_i - l^v_i| > |x_i|, l^v_i é tratado como se fosse infinito, e, da mesma forma, se |x_i - u^v_i| > |x_i|, u^v_i é tratado como se fosse infinito.

Observe que handle_some_primal_gradients_on_finite_bounds_as_residuals não afeta as iterações calculadas, mas afeta apenas o objetivo duplo e os resíduos usados em testes de encerramento e estatísticas informadas.

Redimensionamento

Suponha que recebemos uma matriz de redimensionamento diagonal (coluna) C e uma matriz de redimensionamento diagonal (linha) R com entradas positivas na diagonal. Ao aplicar o redimensionamento como em ShardedQuadraticProgram::RescaleQuadraticProgram, chegamos ao seguinte problema transformado:

\begin{align} \min_{\tilde x} & \, (Cc)^T{\tilde x} + \frac{1}{2}{\tilde x}^T(CQC){\tilde x} \\ \text{s.t.}\; & Rl^{c} \le RAC\tilde x \le Ru^{c} \\ & C^{-1}l^{v} \le \tilde x \le C^{-1}u^{v} \end{align}

Uma solução para o problema original é recuperada como x = C\tilde x. Se \tilde x é uma solução dupla e \tilde r são custos reduzidos para o problema transformado, y = R\tilde y é uma solução dupla e r = C^{-1}\tilde r são reduzidos custos do problema original (derivação omitida).

Identificação da inviabilidade

Um certificado de inviabilidade primaria é um ponto (y, r) \in \mathbb{R}^m \times \mathbb{R}^n satisfatório:

\begin{equation} \left((l^{c})^T[y]_+ - (u^{c})^T[y]_- \right) + \left((l^{v})^T[r]_+ - (u^{v})^T[r]_- \right) > 0 \\ -A^T y = r . \end{equation}

A existência de tal ponto implica que o problema primário não tem uma solução.

Da mesma forma, um certificado de inviabilidade dupla é um ponto x \in \mathbb{R}^n que:

\begin{equation} Qx = 0 \\ c^T x < 0 \\ (Ax)_i \begin{cases} = 0 & \text{if }l^c_i , u^c_i \in \mathbb{R}, \\ \geq 0 & \text{if }l^c_i \in \mathbb{R}, u^c_i = \infty, \\ \leq 0 & \text{if }l^c_i = -\infty, u^c_i \in \mathbb{R}, \end{cases} \\ x_i \begin{cases} = 0 & \text{if }l^v_i , u^v_i \in \mathbb{R}, \\ \geq 0 & \text{if }l^v_i \in \mathbb{R}, u^v_i = \infty, \\ \leq 0 & \text{if }l^v_i = -\infty, u^v_i \in \mathbb{R}. \end{cases} \end{equation}

Para receber certificados de um programa linear, defina Q=0.


  1. Uma matriz de objetivo semidefinita e simétrica S pode ser convertida nesse formato ao fatorar S como S = R^T R (por exemplo, uma fatoração de Cholesky), introduzindo variáveis extras z definidas pelas restrições R x - z = 0, para que x^T S x = z^T z