Fast iterative approximate order-oblivious Orthogonalization algorithm?

Question

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.

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.

Answer 1

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$$

Then

$$ 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}$$

Answer 2

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

Answer 3

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.