Finding Toeplitz matrix nearest to a given matrix For an arbitrary $N\times N$ Hermitian matrix $A$, I want to derive a Toeplitz matrix from $A$ such that the eigenvectors of both matrices have minimal change. 
Specifically I want find the Toeplitz matrix such that the $L^2$ norm between the eigenvectors of the Toeplitz matrix and eigenvectors of the matrix $A$ is minimal. Is there any alternative method other than searching numerically for the matrix? What is the computational cost of such such search?
I am aware of some work done related to perturbations of Toeplitz matrices, in addition eigenvectors of banded toeplitz matrix is studied, but the matrix I want in my application is not banded. I would appreciate any suggestion.
Edit: Is the problem tractable/solvable/realistic if we are given a sequence of matrices $A^n$ instead of $A$?
 A: The set of $n \times n$ symmetric Toeplitz matrices is
$$\left\{ x_1 \mathrm M_1 + x_2 \mathrm M_2 + \cdots + x_n \mathrm M_n  \mid x_1, x_2, \dots, x_n \in \mathbb R \right\}$$
where $\mathrm M_1, \mathrm M_2, \dots, \mathrm M_n$ are $n \times n$ symmetric Toeplitz basis matrices. Let $\mathrm M_1 = \mathrm I_n$ correspond to the main diagonal, whereas the remaining basis matrices correspond to super and sub diagonals.
Let $\mathrm M : \mathbb R^n \to \mbox{Sym}_n (\mathbb R)$ be defined by
$$\mathrm M (\mathrm x) := x_1 \mathrm M_1 + x_2 \mathrm M_2 + \cdots + x_n \mathrm M_n$$

Spectral norm
To complement Suvrit's comment, using the spectral norm, we obtain the following unconstrained optimization problem in $\mathrm x \in \mathbb R^n$
$$\begin{array}{ll} \text{minimize} & \| \mathrm M (\mathrm x) - \mathrm A \|_2\end{array}$$
which can be rewritten as the following semidefinite program (SDP) in $\mathrm x \in \mathbb R^n$ and $t \geq 0$
$$\boxed{\begin{array}{ll} \text{minimize} & t\\ \text{subject to} & - t \,\mathrm I_n \preceq - \mathrm A + x_1 \mathrm M_1 + x_2 \mathrm M_2 + \cdots + x_n \mathrm M_n \preceq t \,\mathrm I_n\end{array}}$$
whose solution can be found numerically.

Frobenius norm
To complement Federico's comment, using the squared Frobenius norm, we obtain the following unconstrained quadratic program (QP) in $\mathrm x \in \mathbb R^n$
$$\begin{array}{ll} \text{minimize} & \| \mathrm M (\mathrm x) - \mathrm A \|_{\text{F}}^2\end{array}$$
where the objective function is
$$\begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_n\end{bmatrix}^\top \begin{bmatrix} \langle \mathrm M_1, \mathrm M_1 \rangle & \langle \mathrm M_1, \mathrm M_2 \rangle & \cdots & \langle \mathrm M_1, \mathrm M_n \rangle\\ \langle \mathrm M_2, \mathrm M_1 \rangle & \langle \mathrm M_2, \mathrm M_2 \rangle & \cdots & \langle \mathrm M_2, \mathrm M_n \rangle\\ \vdots & \vdots & \ddots & \vdots\\ \langle \mathrm M_n, \mathrm M_1 \rangle & \langle \mathrm M_n, \mathrm M_2 \rangle & \cdots & \langle \mathrm M_n, \mathrm M_n \rangle\\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_n\end{bmatrix} - 2 \,\begin{bmatrix} \langle \mathrm A, \mathrm M_1 \rangle\\ \langle \mathrm A, \mathrm M_2 \rangle\\ \vdots \\ \langle \mathrm A, \mathrm M_n \rangle\end{bmatrix}^\top \begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_n\end{bmatrix} + \| \mathrm A \|_{\text{F}}^2$$
where $\langle \mathrm M_i, \mathrm M_j \rangle$ denotes the Frobenius inner product of (symmetric) basis matrices $\mathrm M_i$ and $\mathrm M_j$.
Computing the gradient of the objective function and finding where it does vanish, we obtain the following linear system
$$\begin{bmatrix} \langle \mathrm M_1, \mathrm M_1 \rangle & \langle \mathrm M_1, \mathrm M_2 \rangle & \cdots & \langle \mathrm M_1, \mathrm M_n \rangle\\ \langle \mathrm M_2, \mathrm M_1 \rangle & \langle \mathrm M_2, \mathrm M_2 \rangle & \cdots & \langle \mathrm M_2, \mathrm M_n \rangle\\ \vdots & \vdots & \ddots & \vdots\\ \langle \mathrm M_n, \mathrm M_1 \rangle & \langle \mathrm M_n, \mathrm M_2 \rangle & \cdots & \langle \mathrm M_n, \mathrm M_n \rangle\\ \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_n\end{bmatrix} = \,\begin{bmatrix} \langle \mathrm A, \mathrm M_1 \rangle\\ \langle \mathrm A, \mathrm M_2 \rangle\\ \vdots \\ \langle \mathrm A, \mathrm M_n \rangle\end{bmatrix}$$
Fortunately, the basis matrices are orthogonal and, thus, the matrix above is diagonal. Hence,
$$\begin{bmatrix} n & 0 & \cdots & 0\\ 0 & 2(n-1) & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 2 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_n\end{bmatrix} = \begin{bmatrix} \sum_{i=1}^n a_{i,i}\\ 2 \sum_{i=1}^{n-1} a_{i,i+1}\\ \vdots \\ 2 a_{n,n}\end{bmatrix}$$
and, thus, the solutions $x_1, x_2, \dots, x_n$ are the arithmetic means of the $n$ (distinct) diagonals of $\rm A$
$$\begin{array}{rl} x_1 &= \dfrac{\langle \mathrm A, \mathrm M_1 \rangle}{\langle \mathrm M_1, \mathrm M_1 \rangle} = \dfrac 1n \displaystyle\sum_{i=1}^n a_{i,i}\\ x_2 &= \dfrac{\langle \mathrm A, \mathrm M_2 \rangle}{\langle \mathrm M_2, \mathrm M_2 \rangle} = \dfrac 1{n-1} \displaystyle\sum_{i=1}^{n-1} a_{i,i+1}\\ x_3 &= \dfrac{\langle \mathrm A, \mathrm M_3 \rangle}{\langle \mathrm M_3, \mathrm M_3 \rangle} = \dfrac 1{n-2} \displaystyle\sum_{i=1}^{n-2} a_{i,i+2}\\ & \qquad\qquad\quad\vdots\\ x_n &= \dfrac{\langle \mathrm A, \mathrm M_n \rangle}{\langle \mathrm M_n, \mathrm M_n \rangle} = \displaystyle\sum_{i=1}^1 a_{i,i+n-1} = a_{1,n}\end{array}$$
Lastly, the Toeplitz matrix nearest to the given symmetric matrix $\rm A$ is
$$\boxed{\hat{\mathrm X} := \left( \dfrac 1n \displaystyle\sum_{i=1}^n a_{i,i} \right) \mathrm I_n + \left( \dfrac 1{n-1} \displaystyle\sum_{i=1}^{n-1} a_{i,i+1} \right) \mathrm M_2 + \left( \dfrac 1{n-2} \displaystyle\sum_{i=1}^{n-2} a_{i,i+2} \right) \mathrm M_3 + \cdots + a_{1,n} \mathrm M_n}$$
