I have set of N m-dimensional vectors $\{\phi_i\}$ which gradually loses mutual orthogonality in an algorithm. => I have to re-orthogonalize them every few iterations. But if I use the Gram–Schmidt algorithm for instance, the first vectors (pivot) are perfectly preserved while the the last vectors are strongly perturbed (losing information about their orientation). i.e. the result is sensitive to the order in which I pick the vectors.

I need something which tries to preserve all vectors as much as possible. Like if each has its own inertia.

I was rather thinking to calculate something like orthogonalization forces:

$\vec f_i = -k \sum_{\forall j \neq i} <\vec \phi_i|\vec \phi_j> \vec \phi_j$

($k$ being some stiffness constant ) and run some minimization e.g. gradient descent ($ \phi_i \leftarrow \phi_i+ \vec f_i dt$) or dynamical quenching, possibly even BFGS.

However this is perhaps quite inefficient as each iteration needs $O(mN^2)$ operations, and convergence is quite slow.

Is there more efficient way?

(either using less ops per iteration or with faster convergence)

... but with this property that it preserve orientation of all vectors equally well.


to illustrate the problem for N=3:

  • start with set of vectors $\{ (1,0,0),(0,1,0),(0,0,1) \}$
  • the vectors get perturbed by some process resulting in $\{(1,e_1,e_2),(e_3,1,e_4),(e_5,e_6,1) \}$, where small random numbers $e$ is error caused by the perturbation.
  • Grand-Schmidt will keep direction of $\phi_1 = (\sqrt{1-e_1^2-e_1^2},e_1,e_2)$ and rotate $\phi_2,\phi_3$ acordingly. This means that the whole coordinate system get rotated by $\sin^{-1}(e_1),\sin^{-1}(e_2)$.
  • This can introduces bias or drift, it does not preserve angular momentum. This is a huge problem in many physical applications.
  • For example if $\phi_i$ represent exapnsion coeficients of bonding molecular orbital in atomic orbital basiset (LCAO) e.g. ($s,p_x,p_y,p_z$) this drift introduce fake forces which makes the whole molecule rotate without real physical reason.
  • $\begingroup$ Can you provide context about the application where your system appears? $\endgroup$
    – nicoguaro
    Commented Aug 14, 2019 at 22:21
  • $\begingroup$ kind of Molecular dynamics hybrid between quantum/classical ... $\phi$ represent bonding orbitals between atoms. The drift makes molecule to rotate. $\endgroup$ Commented Aug 15, 2019 at 8:15
  • $\begingroup$ How do your vectors get perturbed? What is commonly done in the literature on your field? $\endgroup$
    – nicoguaro
    Commented Aug 15, 2019 at 10:03
  • 3
    $\begingroup$ This is from 30 year old clouded memories so a comment rather than an answer. Is using SHAKE a la Car-Parrinello MD of any use? If I remember correctly this factors the matrix into the product of a unitary matrix and a symmetric one, so is different from the QR you are using at the moment. But gut feeling is that this is going to be quite tricky to solve ... You might also look at rigid body dynamics in classical MD. $\endgroup$
    – Ian Bush
    Commented Aug 15, 2019 at 11:07
  • $\begingroup$ yes, SHAKE a la Car-Parrinello is one solution for this problem, although not sure how efficient. I was thinking what by formulating the question in more general terms I get bigger picture. $\endgroup$ Commented Aug 15, 2019 at 11:53

3 Answers 3


As already pointed out by Nick Alger, your problem can be restated in a matrix terms. Let's say that $M$ is a matrix where columns are your vectors $\phi_i$. It's possible to prove a theorem that if $m \leq N$ then solution of your problem are columns of matrix

$$Q = M\cdot(M^T\cdot M)^{-1/2}$$

In your vectors were already almost orthogonal, the matrix $M^T\cdot M\approx I$. Let's introduce residual matrix

$$Y = M^T\cdot M - I$$


$$ Q = M \cdot (I+Y)^{-1/2}$$

This expression can be expanded into Taylor series

$$ Q = M \cdot (I - \frac{1}{2}Y+\frac{3}{8}Y^2-\frac{5}{16}Y^3+...)$$

With $Y$ being small, this series converges quickly. In many cases just a first term-approximation will produce good result:

$$Q_1 = M - \frac{1}{2}M \cdot (M^T\cdot M - I)$$

In terms of original vectors, this matrix expression can be rewritten into

$$\hat\phi_i = \frac{3}{2}\phi_i - \frac{1}{2}\sum\limits_j \phi_j \left<\phi_i, \phi_j \right>$$

And that's it - here's your symmetric Gram-Schmidt process, simple and beautiful. Remember that unlike real Gram-Schmidt this one will work only when $\phi_i$ are already almost orthogonal.

All formulas listed above is designed for orthonormal case where we expect that $||\phi_i||=1$ but same solution can be easily adapted to generic orthogonal case by scaling all vectors to unit length before applying same logic. The orthogonalization method in that case would be

$$\hat\phi_i = \phi_i - \frac{1}{2}\sum\limits_{j\neq i} \phi_j \frac{\left<\phi_i, \phi_j \right>}{||\phi_j||^2}$$

  • $\begingroup$ Then apply a QR decomposition to remove the remaining error in the orthonormality. Correct the signs so that the R factor has a positive diagonal. $\endgroup$ Commented Jun 30, 2021 at 13:10
  • $\begingroup$ No, it's not necessary. The beauty of this method is that is works with just matrix multiplications. Or just with vectors. Forget SVD, forget QR: you don't need them if you already have near-orthogonal matrix. I used it in practical CAD system for managing roundoff errors in 3x3 rotation matrices. For a cost of just 3 matrix multiplications I got a solution that was 10 times more accurate than SVD orthogonalization and 200 times faster $\endgroup$ Commented Jun 30, 2021 at 15:39

Depending on how bad the later vectors get, there is a result from Kahan and Parlett that more or less says that doing Gram-Schmidt once, then doing it again to the output of the first pass will sufficiently reduce the residual to low enough tolerances for computer applications. This works really well if you are able to parallelize your Gram-Schmidt algorithm.

If this does not work for you, you could try Modified Gram-Schmidt, which is essentially the same algorithm as Gram-Schmidt but more numerically stable

  • $\begingroup$ Not sure you if this address my problem. I don't speak about numerical stability. I speak about order dependence. $\endgroup$ Commented Aug 14, 2019 at 21:19
  • 1
    $\begingroup$ Order dependence does not matter in exact arithmetics, so this is all about rounding errors and numerical stability. $\endgroup$ Commented Sep 1, 2019 at 11:47
  • $\begingroup$ It does matter - it does affect which vector stay preserved and which vector change. In Gram-Schmidt the first vector is preserved while as you proceede all other vectros move more and more. If the vectors represent something physical which has sort of intertia, this is unphysical process. $\endgroup$ Commented Jul 23, 2021 at 12:05

Let $\Phi$ be the matrix whose columns are the vectors $\phi_i$. You could compute $$\Phi \left(\Phi^T \Phi\right)^{-1/2},$$ which is orthonormal. This is the matrix you would get if you took the SVD, kept all the left and right singular vectors the same, but set the singular values to $1$. As such it is probably the solution of some closest-matrix optimization problem, though I am not sure about this. Anyways, here there is no bias towards any vector.

This reduces the problem to computing a matrix inverse square root. You can do this directly, or iteratively with rational methods. For example, you could combine method 3 in the following paper:

Hale, Nicholas, Nicholas J. Higham, and Lloyd N. Trefethen. "Computing A^α,\log(A), and related matrix functions by contour integrals." SIAM Journal on Numerical Analysis 46.5 (2008): 2505-2523.

with the multishift CG algorithm described in the following paper:

van den Eshof, Jasper, and Gerard LG Sleijpen. "Accurate conjugate gradient methods for families of shifted systems." Applied Numerical Mathematics 49.1 (2004): 17-37.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.