# Fast iterative approximate order-oblivious Orthogonalization algorithm?

I have set of N m-dimensional vectors $$\{\phi_i\}$$ which gradually loose mutual orthogonality in an algorithm. => I have to re-orthogonalize them every few iterations. But if I do e.g. Gram–Schmidt the first vectors (pivot) is perfectly preserved while the the last vectors are strongly perturbed (loosing information about their orientation). i.e. the result depends extremly on 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.

EDIT:

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.
• Can you provide context about the application where your system appears? – nicoguaro Aug 14 at 22:21
• kind of Molecular dynamics hybrid between quantum/classical ... $\phi$ represent bonding orbitals between atoms. The drift makes molecule to rotate. – Prokop Hapala Aug 15 at 8:15
• How do your vectors get perturbed? What is commonly done in the literature on your field? – nicoguaro Aug 15 at 10:03
• 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. – Ian Bush Aug 15 at 11:07
• 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. – Prokop Hapala Aug 15 at 11:53

## 2 Answers

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

• Not sure you if this address my problem. I don't speak about numerical stability. I speak about order dependence. – Prokop Hapala Aug 14 at 21:19
• Order dependence does not matter in exact arithmetics, so this is all about rounding errors and numerical stability. – Jakub Klinkovský Sep 1 at 11:47

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.