Finding the nearest matrix with real eigenvalues

I have no idea what is going on, but your conjecture is not correct. This is more transparent perhaps in the complex version. Consider $$ A = \begin{pmatrix} i & a \\ 0&-i \end{pmatrix} , \quad a>0. $$ This matrix is already in (complex) Schur form, and the obvious procedure to make the eigenvalues real and similar in spirit to what you propose would be to make the diagonal entries real (that is, set them equal to zero here). However, this is not optimal: we can also achieve real spectrum by setting $a_{21}=b$, $b>1/a$, and in fact this is a very small perturbation if $a$ was large. (So the general message is that even eigenvalues that aren't close to the real axis might need only a small perturbation to get them there.)

A real version of this example also works, it just gets somewhat messier computationally. Take $$ A = \begin{pmatrix} 0 & 1 & 0 & 0 \\ -1 & 0 & a & 0\\ 0&0&0&1 \\ 0&0&-1&0 \end{pmatrix} ,\quad a>0 . $$ Your proposal is to remove the $-1$'s, but we can also get real eigenvalues by putting $a_{32}=b$ as long as $ab>4$.

Another counter example to the question which shows how unstable complex eigenvalues can behave: Let $$A=\left(\begin{matrix}0&1&0&\dots &0\\ 0 &0 & 1 &\dots & 0 \\ \vdots &\ddots &\ddots &\ddots& \vdots\\ 0 & \dots&\dots &0 &1\\ 2^{-n} &0 &\dots & \dots &0 \end{matrix}\right)$$ be the $n\times n$ matrix with ones in the first superdiagonal, $2^{-n}$ in the lower left entry and zeros otherwise. Then for large $n$ the matrix $A$ is exponentially close to a nilpotent matrix $B$ (with the lower left entry set to zero) which therefore has only zero as an eigenvalue, in particular purely real spectrum. The eigenvalues of $A$, however, are uniformly distributed on the unit circle. This means that an exponentially small modification of the matrix results in a huge jump in the spectrum. The conjectured closest matrix $A_s$ though, is in spectral norm $\lVert A_s-B\rVert\sim 1$ far away from $B$.

Vanni Noferini and I have published an arxiv e-print, Nearest $\Omega$-stable matrix via Riemannian optimization, in which we address also this problem (in Section 7.5).

Quick summary of our findings:

• There does not seem to be a simple closed-form solution. In particular, you can't get the solution by just truncating the Schur form or the eigendecomposition of $A$.
• We suggest an interesting reformulation of the problem. Namely, let $\operatorname{triu}$ be the operator that returns the upper triangular part of a matrix (diagonal included) $$ \operatorname{triu}(A)_{ij} = \begin{cases} A_{ij} & i\leq j,\\ 0 & i > j \end{cases} $$ and $\operatorname{tril}(A) = A - \operatorname{triu}(A)$ the one that returns its lower triangular part (diagonal excluded). Then the nearest matrix to $A$ with real eigenvalues (returned already in Schur decomposition) is $B = Q\operatorname{triu}(Q^*AQ) Q^*$, where the orthogonal matrix $Q$ solves $$ Q = \arg\min_{Q \in O_n} \|\operatorname{tril}(Q^*AQ)\|_F^2. $$ In this reformulation, the objective function is a nice simple smooth function: it just the sum of squares of the strictly lower triangular entries of $Q^*AQ$. In addition, this formulation has a smaller dimension (the manifold of orthogonal matrices has dimension $\frac{n(n-1)}{2}$ rather than $n^2$) and simpler constraints ("$Q$ is orthogonal" is easier to deal with numerically than "all eigenvalues of $B$ are real").
• Unfortunately this problem (in both formulations, the original one and the alternative formulation) is non-convex and has multiple local minimizers. In similar problems, the number of minimizers grows exponentially with the dimension $n$, and we believe that this holds also for this problem.
• However, code for optimization on matrix manifolds can be used to compute a (local) minimizer numerically in a very fast way, at least for small matrices (up to $n\approx 50$). Our Matlab code is available.
• None of the minimizers suggested by the previous answers for the $3\times 3$ and $4\times 4$ examples in this thread are global minima; we can beat them all. For matrices that are so small, repeating the minimization procedure with various starting points one gets only few local minima, so numerical evidence suggests that the matrices that we compute in Section 7.5 are global minimizers. For instance, for the matrix $$ A = \begin{bmatrix} 1 & 1 & 0\\ -1 & 0 & 0\\ 0 & 0 & 0 \end{bmatrix} $$ in the original question, our method computes the minimizer $$ B \approx \begin{bmatrix} 1.2110 & 0.7722 & -0.0962\\ -0.7722 & -0.2416 & -0.0962\\ -0.0962 & 0.0962 & 0.0306\\ \end{bmatrix}, $$ with $\|A-B\|_F = 0.4946$.
• It is not hard to prove that for each (local) minimizer $B$, $\operatorname{Tr} B = \operatorname{Tr} A$; otherwise one could construct a closer matrix $B + (\operatorname{Tr} A - \operatorname{Tr} B)\frac{1}{n}I$.
• Matrices with real distinct eigenvalues have a neighbourhood formed of matrices with real distinct eigenvalues: hence a (local) minimizer $B$ cannot have distinct eigenvalues. Numerically, often the minimizers have eigenvalues with high multiplicities. For instance, in the example above $B$ has a triple eigenvalue in $1/3$.
• It does not seem easy to extend the results to other matrix norms.
• We started our research by thinking about this problem, but then we were happily surprised to find out that the technique can be applied also to a harder problem that had already been studied in the literature on numerical linear algebra and control theory, that of finding the nearest Hurwitz stable matrix. And, more generally, it can be extended to solve numerically the problem of finding $$ B = \arg \min_{S_\Omega} \|B-A\|_F, $$ where $S_\Omega$ is the set of all (real or complex) matrices with all eigenvalues in a given closed set $\Omega$. Another nice example of how procrastinating and thinking about Mathoverflow questions can lead to useful research sometimes. :)