NOTES | HOME
$$ \newcommand{\RR}{\mathbb{R}} \newcommand{\GG}{\mathbb{G}} \newcommand{\PP}{\mathbb{P}} \newcommand{\PS}{\mathcal{P}} \newcommand{\SS}{\mathbb{S}} \newcommand{\NN}{\mathbb{N}} \newcommand{\ZZ}{\mathbb{Z}} \newcommand{\CC}{\mathbb{C}} \newcommand{\HH}{\mathbb{H}} \newcommand{\ones}{\mathbb{1\hspace{-0.4em}1}} \newcommand{\alg}[1]{\mathfrak{#1}} \newcommand{\mat}[1]{ \begin{pmatrix} #1 \end{pmatrix} } \renewcommand{\bar}{\overline} \renewcommand{\hat}{\widehat} \renewcommand{\tilde}{\widetilde} \newcommand{\inv}[1]{ {#1}^{-1} } \newcommand{\eqdef}{\overset{\text{def}}=} \newcommand{\block}[1]{\left(#1\right)} \newcommand{\set}[1]{\left\{#1\right\}} \newcommand{\abs}[1]{\left|#1\right|} \newcommand{\trace}[1]{\mathrm{tr}\block{#1}} \newcommand{\norm}[1]{ \left\| #1 \right\| } \newcommand{\argmin}[1]{ \underset{#1}{\mathrm{argmin}} } \newcommand{\argmax}[1]{ \underset{#1}{\mathrm{argmax}} } \newcommand{\st}{\ \mathrm{s.t.}\ } \newcommand{\sign}[1]{\mathrm{sign}\block{#1}} \newcommand{\half}{\frac{1}{2}} \newcommand{\inner}[1]{\langle #1 \rangle} \newcommand{\dd}{\mathrm{d}} \newcommand{\ddd}[2]{\frac{\partial #1}{\partial #2} } \newcommand{\db}{\dd^b} \newcommand{\ds}{\dd^s} \newcommand{\dL}{\dd_L} \newcommand{\dR}{\dd_R} \newcommand{\Ad}{\mathrm{Ad}} \newcommand{\ad}{\mathrm{ad}} \newcommand{\LL}{\mathcal{L}} \newcommand{\Krylov}{\mathcal{K}} \newcommand{\Span}[1]{\mathrm{Span}\block{#1}} \newcommand{\diag}{\mathrm{diag}} \newcommand{\tr}{\mathrm{tr}} \newcommand{\sinc}{\mathrm{sinc}} \newcommand{\cat}[1]{\mathcal{#1}} \newcommand{\Ob}[1]{\mathrm{Ob}\block{\cat{#1}}} \newcommand{\Hom}[1]{\mathrm{Hom}\block{\cat{#1}}} \newcommand{\op}[1]{\cat{#1}^{op}} \newcommand{\hom}[2]{\cat{#1}\block{#2}} \newcommand{\id}{\mathrm{id}} \newcommand{\Set}{\mathbb{Set}} \newcommand{\Cat}{\mathbb{Cat}} \newcommand{\Hask}{\mathbb{Hask}} \newcommand{\lim}{\mathrm{lim}\ } \newcommand{\funcat}[1]{\left[\cat{#1}\right]} \newcommand{\natsq}[6]{ \begin{matrix} & #2\block{#4} & \overset{#2\block{#6}}\longrightarrow & #2\block{#5} & \\ {#1}_{#4} \hspace{-1.5em} &\downarrow & & \downarrow & \hspace{-1.5em} {#1}_{#5}\\ & #3\block{#4} & \underset{#3\block{#6}}\longrightarrow & #3\block{#5} & \\ \end{matrix} } \newcommand{\comtri}[6]{ \begin{matrix} #1 & \overset{#4}\longrightarrow & #2 & \\ #6 \hspace{-1em} & \searrow & \downarrow & \hspace{-1em} #5 \\ & & #3 & \end{matrix} } \newcommand{\natism}[6]{ \begin{matrix} & #2\block{#4} & \overset{#2\block{#6}}\longrightarrow & #2\block{#5} & \\ {#1}_{#4} \hspace{-1.5em} &\downarrow \uparrow & & \downarrow \uparrow & \hspace{-1.5em} {#1}_{#5}\\ & #3\block{#4} & \underset{#3\block{#6}}\longrightarrow & #3\block{#5} & \\ \end{matrix} } \newcommand{\cone}[1]{\mathcal{#1}} $$

Linear Complementarity Problem

A few notes on Linear Complementarity (LCP) algorithms for Symmetric Positive Definite (SPD) matrices.

Problem

Given an \(n\times n\) SPD matrix \(M \succ 0\) and a vector \(q \in \RR^n\), the goal is to find \(x \in \RR^n\) such that:

\[0 \leq x \ \bot\ \underbrace{M x + q}_w \geq 0\]

These are the KKT conditions of the following convex QP:

\[\min_{x \geq 0}\ \half x^TMx + q^T x\]

which may be a more suitable point of view for e.g. iterative solvers. Assuming we’re given a Cholesky decomposition of \(M = LL^T\), and letting \(y = L^Tx\) one can rewrite the above as:

\[\min_{L^{-T}y \geq 0} \norm{y - r}^2\]

where \(r = -L^{-1}q\). The KKT conditions become:

\[\exists \lambda \geq 0: \quad y = r + L^{-1} \lambda\]

where \(0 \leq x \ \bot\ \lambda \geq 0\). One may rewrite the cone condition \(L^{-T}y \geq 0\) as \(y \in L^T K\) where \(K = \RR^{n+}\) is the self-dual positive orthant, and letting

\[\mu = L^{-1}\lambda \in L^{-1} K = \block{L^T K}^*\]

we obtain the Moreau decomposition of \(r = y - \mu\) along the cone \(L^T K\) and its (negative) dual. Equivalently, \(q = Mx - \lambda\) is a Moreau decomposition of \(q\) along \(-K\) and its (negative) \(M^{-1}\)-dual.

Projected Jacobi/Gauss-Seidel/SOR

Modulus

Dantzig-Cottle

Lemke

Projected Gradient

Van Bokhoven

Perhaps surprisingly, there exists an explicit formula for solution of the LCP. Of course, there is a catch: the formula uses an exponential number of operations in the size of the problem \(n\).

Let us split the problem using a \(n-1\)-dimensional problem:

\[\mat{M_{11} & M_{1\alpha} \\ M_{\alpha 1} & M_{\alpha \alpha}} \mat{x_1 \\ x_\alpha} + \mat{q_1 \\ q_\alpha} = \mat{\lambda_1 \\ \lambda_\alpha}\]

Now, let us assume we found a solution \(\block{\tilde{x}, \tilde{\lambda}}\) of the \(n-1\)-dimensional LCP \(\block{M_{\alpha\alpha}, q_\alpha}\), we are faced with the following two cases:

We see that in any case, \(\lambda_1 = \block{M_{1\alpha} \tilde{x} + q_1}^+\). This means we can compute \(\lambda_1\) solely from the solution of an \(n-1\)-dimensional problem. This immediately gives rise to an (exponential) algorithm: for each coordinate \(i\), form and solve the corresponding \(n-1\) problem obtained by removing coordinate \(i\), then obtain \(\lambda_i\) from it. While completely untractable for large \(n\) this provides closed-form formula when \(n\) is small, for instance for \(n=2\):

\[\begin{aligned} w_1 &= \block{q_1 + M_{12} \block{-\frac{q_2}{M_{22}}}^{+}}^{+} \\ w_2 &= \block{q_2 + M_{21} \block{-\frac{q_1}{M_{11}}}^{+}}^{+} \\ \end{aligned}\]