# Efficiently removing projection to subspace without having an orthogonal basis

I have a number of vectors $v_1, …, v_n$ and another vector $w$, that all are linearly independent, but not orthogonal. Let $V := \mathop{\text{span}}(v_1, …, v_n)$. I need to remove $w$’s projection to $V$ from $w$, i.e., I need to find the vector $\tilde{w}$ such that $\langle \tilde{w},v \rangle = 0 ~∀~ v∈ V$ and $w-\tilde{w} ∈ V$. I need to do this exactly once for a given $V$.

Note that the scalar product is not the canonical scalar product and has a higher cost than $n$. Touching the scalar product would be extremely tedious, so for the purpose of this question you can consider it a blackbox.

My best solution for this so far is to apply Gram–Schmidt orthonormalisation to the collection $v_1, …, v_n, w$, without normalising the last vector, which then is $\tilde{w}$. With other words, I first orthonormalise $v_1, …, v_n$ to $\hat{v}_1, …, \hat{v}_n$ and then separately remove $w$’s projections to $\hat{v}_1, …, \hat{v}_n$ to obtain $\tilde{w}$.

Doing this, I need to calculate $\frac{n(n-1)}{2}$ scalar products (or norms) just to orthonormalise $V$ as a prerequisite for the last step, which only involves $n$ scalar products. This feels like there could be some more efficient way to do this, but I cannot find one.

• Isn't this just the QR decomposition of $V$? The two other common methods are with Householder reflections and Givens rotations, but in any case you'd probably want to use a standard library implementation of QR. Also, the QR takes $O(n^3)$ time, but the projection takes $O(n^2)$ time, which should make it negligible. – Kirill Nov 30 '16 at 18:48
• Related post on MathOverflow: mathoverflow.net/questions/145688/… – Christian Clason Nov 30 '16 at 18:55
• @Kirill: Yes, orthonormalising $V$ is somewhat equivalent to the $QR$ decomposition. And it’s exactly what I do not need – I only want to project a single vector ($w$). Also, if my scalar product had a cost of $\mathcal{O}(n)$, my $\frac{n(n-1)}{2}$ vs. $n$ would turn into your $\mathcal{O}(n^3)$ vs $\mathcal{O}(n^2)$ – Wrzlprmft Nov 30 '16 at 19:55
• @ChristianClason: It’s a complicated function-based scalar product based on a piecewise Hermite interpolation of a set of unevenly spaced anchor points. Of course, I can write down some spd matrix for it (you can do this for every scalar product on a finite space, IIRC), but that’s already rather tedious. Due to this and some other reasons that need a lot of explanation, you can assume the scalar product to be a blackbox. – Wrzlprmft Dec 1 '16 at 8:45
• So if I understand correctly, your question boils down to "Can I save some effort in (modified) Gram--Schmidt if I use the orthogonalized basis vectors immediately for a projection and don't intend to keep them?" (My guess is no, not in the sense that you can reduce $O(n^3)$ to $O(n^s)$ for $s<3$, but you might be able to lower the constant -- unless there's some special structure in the scalar product and/or the vectors $v$ that you can exploit.) – Christian Clason Dec 1 '16 at 18:23

To fix dimensions and notation, let us say that $V$ is a $m\times n$ matrix with the vectors $v_i$ as columns.

As Jan notes, we have to compute $$\tilde w = [I - V(V^TMV)^{-1}V^TM]w,$$ where $M$ is the symmetric $m\times m$ matrix associated to the scalar product.

If it were possible to solve this problem with $O(n)$ black-box calls to the scalar product $M$, like the OP asks, then it would mean that we can solve the version of the problem in which $M=I$ in $O(mn)$, since now each scalar product costs $O(m)$. This seems too much to ask for, since the usual algorithms to compute projections (QR, SVD, solving the normal equations $(V^TV)^{-1}V^Tw$ explicitly...) all cost $O(mn^2)$.

This answer is an amended and corrected version of the comments: in the comments I claim that computing the scalar product costs $O(m)$ for each matrix $M$, but on second thought this is not true in general. For $M=I$ it's fine, though.

What about simply computing the projection of $w$ along the complement of $V$? $$\tilde w = [I - V(VV^T)^{-1}V^T]w,$$ where $V$ is the matrix that has the basis vectors of $V$ as columns. This $\tilde w$ has all the properties that you look for.

If you are in the non standard scalar product, then there is a symmetric positive matrix $M$ that induces this scalar product. In this case the projection reads $$\tilde w = [I - V(V^TMV)^{-1}V^TM]w.$$

It is OK, if the scalar product is a black box. You only need it's realization.

EDIT: The costs are mainly in the solve with $(V^TMV)$. If CG is applied, there will be a number of scalar product evaluations needed that scales with $n^2$: $\mathcal O(n)$ evaluations times $\mathcal O(n)$ iterations. (read the comments)

• Can you elaborate of the computational complexity of this method? I have the feeling that it would be much worse, given that I have to solve a linear equation system and such. – Wrzlprmft Dec 1 '16 at 10:24
• Depends on how you will do the solve with $(V^TMV)$. In any case it will be less than the orthogonalization of $V$ with respect to $M$. Since $(V^TMV)$ is symmetric positive definite, you can apply CG. With CG, the cost of the overall procedure (and the number of evaluations of the scalar product) scales linearly with $n$. – Jan Dec 1 '16 at 11:19
• Unless I am mistaken, wouldn’t I already need $\frac{n(n+1)}{2}$ evaluations of the scalar product to determine $M$? – Wrzlprmft Dec 1 '16 at 13:03
• Your formula is formally identical to a singular value decomposition and project on the singular vectors (which are an orthonormal set of vectors for the given subspace). Performing the SVD (which is often done via the QR algorithm) is usually numerically more stable than inverting $V^TV$. And the (order of) complexity is for both the same as for Gram-Schmidt. – davidhigh Dec 1 '16 at 13:20
• Are you sure that the number of scalar products scales linearly with $n$ if one uses CG? It seems to me that a linear number of scalar products is required to compute $(V^TMV)x$ for a vector $x$, simply because the output has $n$ entries. So the total number of evaluations is still quadratic. – Federico Poloni Dec 1 '16 at 13:48