\documentclass{slides} \newtheorem{defi}{Definition} \newtheorem{example}{Example} \newtheorem{prop}{Proposition} \newtheorem{lem}{Lemma} \newtheorem{thm}{Theorem} \newtheorem{cor}{Corollary} \newtheorem{rem}{Remark} \newtheorem{conj}{Conjecture} \newtheorem{alg}{Algorithm} \newtheorem{ex}{Exercise} \newcounter{count} \newcommand{\adj}{{\rm adj\,}} \newcommand{\trace}{{\rm trace\,}} \newcommand{\diag}{{\rm diag\,}} \newcommand{\Diag}{{\rm Diag\,}} \newcommand{\conv}{{\rm conv\,}} \newcommand{\bs}{\setcounter{count}{0} \section} %\addtolength{\topmargin}{-.5in} %\addtolength{\textheight}{.5in} %\addtolength{\oddsidemargin}{.25in} %\addtolength{\evensidemargin}{.25in} \def\QED{~\rule[-1pt] {8pt}{8pt}\par\medskip ~~} \begin{document} \bibliographystyle{plain} \begin{slide}{} \begin{center} {\bf MAX-MIN EIGENVALUE PROBLEMS, PRIMAL-DUAL INTERIOR POINT ALGORITHMS, and TRUST REGION SUBPROBLEMS} \end{center} Rendl-Vanderbei-Wolkowicz, Optimization Methods and Software, 1995, Vol 5, pp 1-16. \begin{center} Outline \end{center} $\bullet$ two Primal-dual interior point algorithms for maximizing the smallest eigenvalue of a symmetric matrix over diagonal perturbations\\ $\bullet$ One of the algorithms does this transformation through an intermediate transformation to a trust region subproblem. This allows the removal of a dense row and results in a nonlinear SDP. \end{slide} \begin{slide}{} {\em max-min eigenvalue problem} \[ ~~~(MMP)~~~\omega^*:=\max_{v^te=0} \lambda_{\min} (C^b - \Diag(v)), \] where $C^b$ is an $(n+1) \times (n+1)$ real symmetric matrix, $v \in \Re^{n+1}$, and $\Diag(v)$ denotes the diagonal matrix formed from the vector $v$. \end{slide} \begin{slide}{} $\bullet$ Applications e.g.: provides bounds in branch and bound codes for graph bisection; many engineering applications. $\bullet$ In this paper we study two primal-dual interior point algorithms which show that these problems can be solved robustly and quickly for very large dimensions. $\bullet$ difficulties: loss of differentiability due to multiplicity of the smallest eigenvalue at the optimum $\bullet$ smallest eigenvalue is a concave function, so subgradient approaches can be used $\bullet$ Newton-based algorithms with local quadratic convergence are still possible though the local convergence depends on correctly identifying the multiplicity of the smallest eigenvalue $\bullet$ just as in ordinary linear programming, the interior point methods we present do not seem to bothered by degeneracy. \end{slide} \begin{slide}{} \begin{center} {\bf NOTATION} \end{center} {\bf real symmetric matrices} $\cal S$\\ {\bf trace inner product} $\left< A,B \right> := \trace AB.$ \\ {\bf positive semidefinite matrices} $\cal P$ or ${\cal S}_+$ forms a closed convex cone, which is {\bf self-polar}, \[{\cal P}^+ := \{ K=K^t: \trace KP \geq 0,~ \forall P \in {\cal P} \} = {\cal P}.\] Lo\"{e}wner partial order, i.e. $A \succ$ (resp. $\succeq$) $B$ denotes $A-B$ is positive definite (resp. semidefinite). $v \in \Re^n,$ matrix $M$:\\ $\Diag(v)$ - diagonal matrix, $\diag(M)$ - column vector\\ $e \in \Re^n$ - {\bf vector of ones}; $e_i$ {\bf i-th unit vector}; $E_i$ matrix with 1 in $i,i$ position\\ $M^\dagger$ - {\bf Moore-Penrose generalized inverse}\\ ${\cal R}(M),~{\cal N}(M)$ denote {\bf range space} and {\bf null space}, respectively\\ $\det(M)$ denotes the {\bf determinant}\\ $M \circ N$ {\bf Hadamard product} or entrywise product, \end{slide} \begin{slide}{} \begin{center} {\bf EQUIVALENT FORMULATIONS} \end{center} subtract elements summing to 0 from the diagonal \[ C^b = \left[ \begin{array}{cc} 0 & -b^t \\ -b & C \end{array} \right]. \] (or choose dense row and column) an equivalent max-min problem is the {\em parametric trust region subproblem} \[ \begin{array}{rclr} ~(TRP)~~~(n+1)\omega^*&=&\max_{~~u^te=0~~} \min_{x^tx=n} \\ && ~~~x^t(C-\Diag(u))x-2b^tx, \end{array} \] where $x,u \in \Re^n.$ \end{slide} \begin{slide}{} from duality an equivalent problem is the (dualTRP) \[ \begin{array}{rclr} (n+1)\omega^*&=&\max_{~~u^te=0~~} \max_{C-\Diag(u)-\lambda I \succeq 0} \\ && ~~~n \lambda - b^t (C-\Diag(u)-\lambda I)^{\dagger} b, \end{array} \] with $y=u+\lambda e$, this is equivalent to the {\bf Nonlinear SDP} \[ (D)~~~~\begin{array}{llll} (n+1) \omega^*= \max & f(y) := \\ &~~~y^te- b^t (C-\Diag(y))^{\dagger} b& \\ {\rm subject~to} & C-\Diag(y) \succeq 0\\ & ~y \in \Re^n, \end{array} \] where $e$ is the vector of ones. \end{slide} \begin{slide}{} {\bf THEOREM A}\\ Suppose that $y$ solves (D) with value $f(y)=(n+1)\omega^*$. Let: \[ \lambda := \frac{y^te}{n},~~u:=y-\lambda e;\] \[ \bar{x} = C_y^{\dagger}b,~~\alpha:=\bar{x}^t\bar{x};\] \[ x=\bar{x}+w,~ \mbox{ for some}~ w \in {\cal N}(C_u) ~\mbox{such that}~ x^tx=n; \] and \[\bar{z}:=\left( \begin{array}{c} 1 \\ \bar{x} \end{array} \right),~~ z:=\left( \begin{array}{c} 1 \\ x \end{array} \right) \] \[ t:=\lambda+b^t\bar{x},~~ v:= \left(\begin{array}{c} \frac n{n+1}t \\ u-\frac t{n+1} e \end{array} \right). \] \end{slide} \begin{slide}{} Then: \begin{enumerate} \item $u,\lambda$ solves (dualTRP) with optimal value $(n+1)\omega^*$; \item $\alpha \leq n$, $w$ exists for above and $u,x$ solves the (TRP) with optimal value $(n+1)\omega^*$, \item $z,\bar{z}$ are eigenvectors for the optimal eigenvalue $\omega^*$ with optimal perturbation $v$ for (MMP). \end{enumerate} \QED \end{slide} \begin{slide}{} a dual problem to problem (D) is the min-max of the Lagrangian \[ \min_{X \succeq 0} \max_y f(y) + \trace X(C-\Diag(y)). \] adding the stationary point for the concave function gives the primal problem, or the Wolfe dual of (D), \[ ~~~(P)~~~~\begin{array}{ll} \omega^*= \min & f(y) + \trace X (C-\Diag(y)) \\ {\rm subject~to} & \nabla f(y) -\diag(X) = 0\\ & X \succeq 0. \end{array} \] no duality gap between primal and dual problems for $y,X$ a feasible pair, duality gap is just complementary slackness \[ \trace X (C-\Diag(y)). \] dual log-barrier with log-barrier parameter $\mu \downarrow 0$: \[ \begin{array}{llll} \max & f(y) +\mu \log \det (C-\Diag(y)) & \\ {\rm subject~to} & ~y \in \Re^n; \end{array} \] \end{slide} \begin{slide}{} \begin{center} {\bf ALTERNATIVE FORMULATION} \end{center} \[ \omega^*=\max \{\omega : v^te=0,~C^b - \Diag(v) \succeq \omega I \}. \] equivalently with $y^s=v+\omega e$ \[ \omega^*=\max \{e^ty^s: ~C^b - \Diag(y^s) \succeq 0 \} \] with dual \[ \begin{array}{llll} \min & \trace C^bX^a& \\ {\rm subject~to} & \diag(X^a)= e \\ & ~X^a\succeq 0. \end{array} \] \end{slide} \begin{slide}{} relationships with previous formulation using the Schur complement\\ {\bf THEOREM B} Let $X^a,y^s$ be an optimal pair for the above dual pair, with optimal value $\omega^*$. Partition \[ X^a := \left[ \begin{array}{cc} 1 & a^t \\ a & W \end{array} \right], ~y^s:=\left[ \begin{array}{c} s \\ y \end{array} \right]. \] Let: \[ v:=y^s-\omega^*e,~t:=\frac{n+1}n (s-\omega^*),~\bar{x}:=a;\] \[ \lambda := t-b^t\bar{x}~~u:=y-\lambda e;\] and $x:=\bar{x}+w,$ for some $w \in {\cal N}(C-\Diag(u))$ such that $x^tx=n.$ Then 1.,2.,3. of Theorem A hold. Moreover, $X:=W-aa^t$ solves (D) the dual program of Theorem A. \QED \end{slide} \begin{slide}{} \begin{center} {\bf DERIVATIVES for Newton's Method} \end{center} \[ f(y) := y^te- b^t (C-\Diag(y))^{-1} b. \] Denote the log-barrier term by \[ g(y) := \log \det (C-\Diag(y)), \] \[ C_y := C - \Diag(y), ~ x_y := C_y^{-1} b. \] From \[ C_y x_y = b, \] we differentiate with respect to $y$ to get \[ \frac{\partial C_y}{\partial y} x_y + C_y \frac{\partial x_y}{\partial y} = 0, \] \end{slide} \begin{slide}{} i.e. the $n \times n$ Jacobian transposed \[ \frac{\partial x_y}{\partial y} = C_y^{-1} \Diag(x_y), \] while \[ \frac{\partial x_y}{\partial y_i} = C_y^{-1} (E_ix_y). \] Therefore \[ \begin{array}{rcl} \frac{\partial f(u)}{\partial y_i} &=& 1 -b^tC_y^{-1} (E_ix_y)\\ &=&1 -||E_ix_y||^2. \end{array} \] This yields the gradient \[ \begin{array}{rcl} \frac{\partial f(u)}{\partial y} &=& e-\diag(x_yx_y^t). \end{array} \] \end{slide} \begin{slide}{} Given a symmetric matrix $B$, \[ \frac{\partial \det B}{\partial B} = \adj(B),\] where $\adj$ denotes the {\em adjoint matrix}. Therefore, by Cramer's rule, \[ \begin{array}{rcl} \frac{\partial g(y)}{\partial y_i} &=& \frac {-1}{\det C_y} \trace E_i \adj(C_y)\\ &=&-\trace E_i C_y^{-1}, \end{array} \] and \[ \begin{array}{rcl} \frac{\partial g(y)}{\partial y} &=& -\diag( C_y^{-1}). \end{array} \] If we differentiate $ C_y x_y = b, $ twice we get \[ -E_i \frac{\partial x_y}{\partial y_j} -E_j \frac{\partial x_y}{\partial y_i} + C_y \frac{\partial^2 x_y}{\partial y_i \partial y_j} =0, \] which yields \[ \begin{array}{rcl} \frac{\partial^2 x_{u}}{\partial y_i \partial y_j} &=& C_y^{-1}( E_i \frac{\partial x_y}{\partial y_j} +E_j \frac{\partial x_y}{\partial y_i}). \end{array} \] \end{slide} \begin{slide}{} Also \[ \begin{array}{rcl} \frac{\partial^2 f(u)}{\partial y_i \partial y_j} &=& -2 (E_i x_y)^t \frac{\partial (E_ix_y)}{\partial y_j}\\ &=&-2 (E_i x_y)^t E_i \frac{\partial x_y}{\partial y_j}\\ &=& -2(E_ix_y)^tC_y^{-1} (E_j x_y). \end{array} \] Therefore, \[ \begin{array}{rcl} \frac{\partial^2 f(u)}{\partial y^2 } &=& -2 x_y x_y^t \circ C_y^{-1}, \end{array} \] where $\circ$ denotes the {\em Hadamard product}. \end{slide} \begin{slide}{} For a nonsingular matrix function $B(t)$, we can differentiate both sides of $B(t)B(t)^{-1} = I$ and obtain the derivative \[ \frac{\partial B(t)^{-1}}{\partial t } = -B(t)^{-1} \frac{\partial B(t)}{\partial t }B(t)^{-1}. \] Therefore, differentiating $g$ twice yields: \[ \begin{array}{rcl} \frac{\partial g(y)}{\partial y_i \partial y_j} &=& -\frac{\partial \,(\trace E_i C_y^{-1})}{\partial y_j}\\ &=&-\trace E_i C_y^{-1}E_j C_y^{-1}, \end{array} \] and \[ \begin{array}{rcl} \frac{\partial g(y)}{\partial y^2 } &=& - C_y^{-1} \circ C_y^{-1}. \end{array} \] \end{slide} \begin{slide}{} In summary, the gradients and Hessians are: \[ \nabla f = e-\diag(x_yx_y^t); \] \[ \nabla g = -\diag( C_y^{-1}); \] \[ \nabla^2 f = -2x_y x_y^t \circ C_y^{-1}; \] \[ \nabla^2 g = -C_y^{-1} \circ C_y^{-1}. \] Note that the above implies that the functions $f,g$ are concave when $C_y$ is positive definite, since Schur's Theorem on Hadamard products states that the cone of positive semidefinite matrices is closed under the Hadamard product. \end{slide} \begin{slide}{} \begin{center} {\bf PRIMAL-DUAL METHODS} \end{center} by adding a slack variable matrix $Z$ \[ ~~~(D)~~~~\begin{array}{llll} \omega^*= \max & f(y) := y^te- b^t (C-\Diag(y))^{\dagger} b& \\ {\rm subject~to} & C-\Diag(y) -Z = 0\\ & ~u \in \Re^n,~Z \succeq 0. \end{array} \] The primal-dual log-barrier problem is then \[ ~~~(B)~~~~\begin{array}{llll} \omega^*_\mu= \max & f(y) + \mu \log \det (Z) &\\ {\rm subject~to} & C-\Diag(y) -Z = 0. \end{array} \] The Lagrangian is \[ \begin{array}{rcl} L(y,Z,X) &=& f(y) + \mu \log \det (Z)\\ &&~~~~ +\trace X(C-\Diag(y) -Z), \end{array} \] diff. yields opt. cond. $F(y,Z,X)=0:$ \[ \begin{array}{rcl} e-\diag(x_yx_y^t) - \diag(X) &=& 0~~\mbox{(primal feasibility)}\\ -X + \mu Z^{-1} &=& 0~~\mbox{(comp. slack.)}\\ C-\Diag(y) -Z &=& 0~~\mbox{(dual feasibility)}\\ Z \succeq 0,~ X \succeq 0. \end{array} \] \end{slide} \begin{slide}{} $\bullet$ Primal feasibility corresponds to feasibility of $x_y$ in (P), i.e. $\diag(X)=0$ corresponds to $x_y$ feasible in TRS.\\ Emphasizes that no feasibility due to the hard case.\\ $\bullet$ However, $ZX=0$ means columns of $X$ and also of $X^{\frac 12}$ form a basis for the null space of $Z=C-\Diag(y)$.\\ $\bullet$ Linear combination of the columns of $X^{\frac 12}$ added to $x_y$ gives a feasible optimal solution to subproblem. \end{slide} \begin{slide}{} apply Newton's method; system $-F^\prime \delta = F$ is: \[ \begin{array}{rcl} \left[ \begin{array}{ccc} 2x_y x_y^t \circ C_y^{-1} & 0 & \diag(\cdot) \\ 0 & \mu Z^{-1} \cdot Z^{-1} & I \\ \Diag(\cdot) & I & 0 \\ \end{array} \right] \left( \begin{array}{c} \delta y \\ \delta Z \\ \delta X \end{array} \right) &=&\\ ~~~~~~ \left( \begin{array}{c} ( e-\diag(x_yx_y^t) - \diag(X))\\ ( -X + \mu Z^{-1})\\ (C-\Diag(y) -Z )\\ \end{array} \right) \end{array} \] first equation \[ \delta Z = Z - \frac 1\mu ZXZ - \frac 1\mu Z(\delta X)Z. \] Moreover, feasibility maintained by linear equations Therefore, \[ \delta Z = - \Diag(\delta y),~C_y^{-1} = Z^{-1}. \] eliminate $\delta Z$ \[ \left[ \begin{array}{cc} -\frac 1\mu Z \cdot Z & \Diag(\cdot) \\ \diag (\cdot) & -\nabla^2 f(y) \end{array} \right] \left( \begin{array}{c} \delta X \\ \delta y \end{array} \right) = \left( \begin{array}{c} -Z+ \frac 1\mu ZXZ \\ \nabla f(y) - \diag(X) \end{array} \right) \] \end{slide} \begin{slide}{} multiplying through by $\mu$ and $Z^{-1}$ and eliminating $\delta X$ \begin{eqnarray*} \left[ \begin{array}{cc} I & -\mu Z^{-1}\Diag(\cdot)Z^{-1} \\ 0 & \diag(\mu Z^{-1}\Diag(\cdot)Z^{-1}) - \nabla^2 f(y) \end{array} \right] \left( \begin{array}{c} \delta X \\ \delta y \end{array} \right) =\\ \left( \begin{array}{c} \mu Z^{-1} - X \\ \diag(-\mu Z^{-1} + X) +\nabla f(y) - \diag(X) \end{array} \right), \end{eqnarray*} so that \[ \delta X = \mu Z^{-1} - X + \mu Z^{-1}\Diag(\delta y)Z^{-1}. \] Summary \begin{eqnarray} \delta y& = &( (2x_y x_y^t + \mu Z^{-1}) \circ Z^{-1} )^{-1} (e- \diag(x_y x_y^t)-\mu \diag(Z^{-1}) );\nonumber \\ \delta Z& =& - \Diag(\delta y); \\ X+\delta X &= &\mu Z^{-1} ( Z + \Diag(\delta y))Z^{-1}\nonumber. \end{eqnarray} check both $X,Z$ positive definite from $ Z \pm \Diag(\delta y)$ \end{slide} \begin{slide}{} \begin{center} {\bf ALTERNATIVE ALGORITHM} \end{center} primal-dual optimality conditions: \[ \begin{array}{rcl} e- \diag(X^a) &=& 0~~\mbox{(primal feasibility)}\\ -X^a + \mu Z^{-1} &=& 0~~\mbox{(compl. slack.)}\\ C^b-\Diag(y^s) -Z &=& 0~~\mbox{(dual feasibility)}\\ Z \succeq 0,~ X^a \succeq 0. \end{array} \] $X^a,Z$ primal and dual matrices but one dimension larger - system like an LP apply Newton's method \[ \begin{array}{rcl} \left[ \begin{array}{ccc} 0 & 0 & \diag(\cdot) \\ 0 & \mu Z^{-1} \cdot Z^{-1} & I \\ \Diag(\cdot) & I & 0 \end{array} \right] \left( \begin{array}{c} \delta y^s \\ \delta Z \\ \delta X^a \end{array} \right) &=&\\ ~~~~~~ \left( \begin{array}{c} ( e- \diag(X^a))\\ ( -X^a + \mu Z^{-1})\\ (C^b-\Diag(y^s) -Z ) \end{array} \right). \end{array} \] similar equations \end{slide} \begin{slide}{} \begin{center} {\bf LINE SEARCHES} \end{center} For algorithm 1, verify that the new $X$ and $Z$ are positive definite. \[ Z - t \Diag(\delta y) \succ 0,~\mbox{for}~Z.\] $t>0$ step length. For $X$, when the error $||\mu Z^{-1} -X||$ is small, then a reasonable heuristic is to check \[ Z + t \Diag(\delta y) \succ 0,~\mbox{for}~X.\] Let $Z=R^tR$ be the Cholesky factorization of $Z$. Then the above is equivalent to finding $t>0$ such that \[ I \pm t R^{-t}(\Diag(\delta y))R^{-1} \succ 0. \] With this change, the problem is equivalent to the step length problem in ordinary linear programming, i.e. suppose that the matrix $\Diag(\delta y)$ is not negative semidefinite (in which case the step length can be infinite and so is chosen to be 1). Then we get that the step for $Z$ is \[ t_Z := \min \{ 1, \frac 1{\lambda_{\max} (R^{-t}(\Diag(\delta y))R^{-1})} \}; \] while, if the matrix $\Diag(\delta y)$ is not positive semidefinite (in which case the step length can be infinite and so is chosen to be 1) the step for $X$ is \[ t_X := \min \{ 1, \frac 1{-\lambda_{\min} (R^{-t}(\Diag(\delta y))R^{-1})} \}. \] By the above we see that we can use any upper bound for the maximum eigenvalue and any lower bound for the minimum eigenvalue. Note that for an $n \times n$ symmetric matrix $K$, \[ \lambda_{\max} (K) \leq \frac{\trace K}{n} + \sqrt{n-1}\sqrt{ \frac {\trace K^2}{n} - \left( \frac{\trace K}{n}\right)^2 }; \] \[ \lambda_{\min} (K) \geq \frac{\trace K}{n} - \sqrt{n-1}\sqrt{ \frac {\trace K^2}{n} - \left( \frac{\trace K}{n}\right)^2 }. \] We can find the traces using $Z^{-1}$, e.g. \[ \begin{array}{rcl} \trace (R^{-t}(\Diag(\delta y))R^{-1}) &=& \trace \Diag(\delta y) Z^{-1} \\ &=& \delta y^t \diag( Z^{-1}). \end{array} \] Alternatively, to find the bounds, we could first shift $R^{-t}(\Diag(\delta y))R^{-1}$ by a multiple of the identity to guarantee that it is positive definite and then calculate its norm. We then shift by a multiple of the identity to guarantee that it is negative definite and find the norm of the negative of the matrix. By shifting back, we get the values for the largest and smallest eigenvalues. Alternatively, we could calculate the largest and smallest eigenvalues of $\Diag(\delta y) Z^{-1}$ directly. Note that if $t_X=1$, then $\delta X$ is a good step for $X$. However, if $t_X<1$, then $X+t_X\delta X$ is not necessarily equal to $ \mu Z^{-1} ( Z + \Diag(t_X \delta y))Z^{-1}.$ This situation almost never arose in practice. But an extra safeguard was added to ensure that $X$ remain positive definite. After finding these maximum steplengths, we then multiply the step length by .90 to guarantee that we do not get too close to the boundary. Because of the nonlinearity in algorithm 1, we then take one steplength for both variables; while algorithm 2 uses different steplengths in the primal and dual variables. From numerical tests, it appears that the step length 1 for $X$ is usually not too large, i.e. it does not usually lose positive definiteness. But, the step length for $Z$ is much too large immediately after $\mu$ is decreased. When $\mu$ is not decreased, the step length 1 seems good for both variables. An efficient line search could be done using a Lanczos type algorithm to calculate the smallest eigenvalues of $C_y - \alpha \delta y,~ C_y + \beta \delta y$. Or inverse iteration can be used so the matrices from the Cholesky factorization do not get inverted. \end{slide} \begin{slide}{} \begin{center} {\bf CONCLUSION} \end{center} We have derived two primal-dual algorithms for maximizing the minimum eigenvalue of a diagonally perturbed symmetric matrix. Both algorithms have been tested extensively. The results for the first algorithm, which is based on the trust region subproblem, were similar. Therefore, if the matrix has a dense row, the first algorithm should be used since the dense row can be eliminated. \end{slide} \end{document}