On computing the distance to stability for matrices using linear dissipative Hamiltonian systems

In this paper, we consider the problem of computing the nearest stable matrix to an unstable one. We propose new algorithms to solve this problem based on a reformulation using linear dissipative Hamiltonian systems: we show that a matrix A is stable if and only if it can be written as A=(J−R)Q, where J=−JT, R⪰0 and Q≻0 (that is, R is positive semidefinite and Q is positive definite). This reformulation results in an equivalent optimization problem with a simple convex feasible set. We propose three strategies to solve the problem in variables (J,R,Q): (i) a block coordinate descent method, (ii) a projected gradient descent method, and (iii) a fast gradient method inspired from smooth convex optimization. These methods require O(n3) operations per iteration, where n is the size of A. We show the effectiveness of the fast gradient method compared to the other approaches and to several state-of-the-art algorithms.