Is there any fast way to solve a quadratic matrix equation in Mathematica approximately?

One method is to use Mathematica's PseudoInverse[] on the matrix equation that involves the eigenvalues of $M$. To arrive at this matrix equation, multiply each side of $M^2 + \alpha M + \beta I \approx 0$ by an eigenvector $v_1$ of $M$, to get

$$ (\lambda_1^2 + \alpha \lambda_1 + \beta) v_1 \approx 0 \implies \lambda_1^2 + \alpha \lambda_1 + \beta\approx 0,$$ where $\lambda_1$ is the corresponding eigenvalue for $v_1$. Repeat this for all three eigenvectors (assuming that $M$ is semi-simple) to arrive at the following approximate matrix equation

$$\begin{pmatrix}\lambda_1 & 1 \\ \lambda_2 & 1 \\ \lambda_3 & 1 \end{pmatrix} \begin{pmatrix} \alpha \\ \beta\end{pmatrix} \approx - \begin{pmatrix} \lambda_1^2 \\ \lambda_2^2 \\ \lambda_3^2 \end{pmatrix}. $$ To solve this use in Mathematica

{α,β} = PseudoInverse[{{λ1,1},{λ2,1},{λ3,1}}].{λ1^2,λ2^2,λ3^2} 

Then all you need are the eigenvalues of $M$, which you can get through

{λ1,λ2,λ3} = Eigenvalues[M] 

For a very large matrix it maybe impractical to find all of $M$'s eigenvalues, in which case you can use the Power Method to obtain, for example, the three largest eigenvalues and then use the above method to find α and β. The more eigenvalues you use the better the approximation.

Hope this helps.

By construction, M is real symmetric and almost certainly positive definite.
The eigenvalues of M are e = (d/d[[1]])^2, where d is the vector of singular values of A.
The eigenvalues of C = M^2 + alpha M + beta are e^2 + alpha e + beta.
The eigenvectors M and C are the same as the left singular vectors of A, say V.
The offdiagonals of V'CV are zero.
The sum of squares of the diagonals of V'CV = the sum of squares of all the elements of C.
But the diagonals of V'CV are just the eigenvalues of C, so we can simplify the problem:
minimize Tr[(e^2 + alpha e + beta)^2].
You already know e, because you constructed it from d, so the solution is immediate:
{alpha = Covariance[e,-e^2]/Variance[e], beta = Mean[-e^2] - alpha*Mean[e]}.
That's the exact least-squares solution.

EDIT - It occurred to me that, although getting the singular values of a 1000 x 1000 matrix takes only a few seconds, if you have to do it many times (e.g., inside a loop) then the total time may be prohibitive, so further analysis may help.

Simulations suggest that if the elements of A are iid Uniform[-1,1] then the expected value of e is close to Reverse[(#/2 + ArcSin[#]/Pi)^2 & [Range@n/n]]. This code

{n = #, y = -(x = (#/2 + ArcSin[#]/Pi)^2 & [Range@n/n])^2; N@{a = (n x.y - (Tr@x)(Tr@y))/(n x.x - (Tr@x)^2), b = (Tr@y - a*Tr@x)/n}}& /@ {10,100,1000} 

gives {alpha,beta} pairs for n = {10, 100, 1000}.

{{10, {-0.939996, 0.104015}} {100, {-0.764133, 0.0667884}} {1000, {-0.743093, 0.0626846}}} 

The pair for n = 1000 should be a good approximation to the actual least-squares coefficients for a sample 1000 x 1000 matrix.

EDIT 2 - I get the same results when the elements of A are normal and logistic as when they're uniform. Has taking n = 1000 put us in asymptoticland? Would someone please check these results?

EDIT 3 - I think what we're seeing is a result of the quarter-circle law that the asymptotic density of the singular values of an n x n matrix whose elements are iid zero-mean random variables with variance 1/n converges to Sqrt[4 - x^2]/Pi, 0 <= x <= 2. The singular values we are dealing with have been normalized -- i.e., divided by the largest one -- so (rounding some statistical corners) their asymptotic density is approximately proportional to Sqrt[1 - x^2], 0 <= x <= 1, and their cumulative distribution is approximately F[x] := 2(x Sqrt[1 - x^2] + ArcSin[x])/Pi. The function I gave in my first edit, say K[p] := p/2 + ArcSin[p]/Pi, is an approximate inverse of F[x]. The actual inverse has no closed form and must be found numerically, say X[p] := Block[{x},x/.FindRoot[F[x]==p,{x,K[p]}]]. For n = 1000, the rms error in approximating the normalized singular values using X[p] is typically much smaller than the rms error using K[p]; e.g., .00095 vs .0091 .