# Orthonormalize a matrix in place

**URL:** <https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887>\
**Category:** Performance\
**Tags:** linearalgebra\
**Created:** [January 5, 2021, 5:00pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887 "2021-01-05T17:00:18Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![matthiasbe](https://avatars.discourse-cdn.com/v4/letter/m/5fc32e/32.png) [@matthiasbe](https://discourse.julialang.org/u/matthiasbe)\
**Post date:** [January 5, 2021, 5:00pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/1 "2021-01-05T17:00:18Z")

</div>

Hi,

For a given matrix M \in \mathbb{R}^{m \times n}, m \geq n , I would like to orthonormalize the columns of M in place.

M \leftarrow \texttt{qr}(M)\texttt{.Q}

I see that the function `LinearAlgebra.qr!` aims at something similar, but I cannot figure out what will be put in the matrix M (will it be R, or Q or something else ?). From experiments it feels that the following is put

- the R factor on the upper triangular part of the matrix,
- the reflector on the rest

It makes sens because this is how it is done with LAPACK.

Is there a way to put the economic (meaning m \times n) Q factor in M without using extra memory ? I should say that the R factor is not needed in my usecase.

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [January 5, 2021, 5:50pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/2 "2021-01-05T17:50:28Z")

</div>

I am not a numerical analyst, but I would guess that if you aren’t concerned with numerical stability or speed, you could just write your own Gram-Schmidt that updates in-place. Maybe that is literally never advisable, but at least in principle something like that would give you what you want.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 5, 2021, 6:02pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/3 "2021-01-05T18:02:30Z")

</div>

> [@matthiasbe](#):
>
> Is there a way to put the economic (meaning m×nm \times n ) Q factor in M without using extra memory ? I should say that the R factor is not needed in my usecase.

The Householder QR algorithm (used by the `qr` function) does not compute or store the Q matrix explicitly — it uses an implicit representation of Q in terms of a sequence of elementary reflectors. (In most applications, you don’t need Q explicitly, only a way to multiply any given vector or matrix by Q or Q^\*.)

You could implement modified Gram–Schmidt yourself, of course, as @cgeoga suggested, but there is some numerical loss of orthogonality that might bite you unless you use it carefully.

What are you trying to do? What do you need Q for?

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [January 5, 2021, 6:22pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/4 "2021-01-05T18:22:07Z")

</div>

If you orthogonalize with classical Gram-Schmidt twice you’ll be fine. The GMRES code in ```SIAMFANLEquations.jl```` does it that way. That code, and the orthogonalizer live in

[https://github.com/ctkelley/SIAMFANLEquations.jl/tree/master/src/Solvers/LinearSolvers](https://github.com/ctkelley/SIAMFANLEquations.jl/tree/master/src/Solvers/LinearSolvers)

Here is a quick and dirty qr code that uses the orthogonalizer and overwrites `A` with Q.  
The default option, `cgs2` is classical Gram-Schmidt twice, which is the best choice on any modern computer.

```julia
function qrctk!(A, orth = "cgs2")
    T = typeof(A[1, 1])
    (m, n) = size(A)
    R = zeros(T, n, n)
    @views R[1, 1] = norm(A[:, 1])
    @views A[:, 1] /= R[1, 1]
    @views for k = 2:n
        hv = vec(R[1:k, k])
        Qkm = view(A, :, 1:k-1)
        vv = vec(A[:, k])
        Orthogonalize!(Qkm, hv, vv, orth)
    end
    return (Q = A, R = R)
end

```

---

<div class="post-metadata">

**Author:** ![matthiasbe](https://avatars.discourse-cdn.com/v4/letter/m/5fc32e/32.png) [@matthiasbe](https://discourse.julialang.org/u/matthiasbe)\
**Post date:** [January 6, 2021, 8:35am UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/5 "2021-01-06T08:35:50Z")

</div>

Thank you for these answers.

@ctkelley @cgeoga Having already tried to implement the Gram-Schmidt algorithm in the past, I fear I won’t be able to match the stability and the performance of the library’s QR method. Thus my question was toward using the existing function of `LinearAlgebra` to achieve my goal, if possible.

@stevengj Here is a usecase

I want to create an orthogonal matrix distributed on different processors, but this matrix can be very large, so it shouldn’t be generated all at once and then split the columns.

I generate the different partitions on the different processors, and then I orthogonalize the partitions against one another.

```julia
function orthonormalize(A::DArray)
  for p in workers()
    # Compute the local orthonormal partition
    @spawnat p A[:l] = Matrix(LinearAlgebra.qr(A[:l]).Q)

    # On processors p+1 .. P, orthogonalize against the basis of p
    for p2 in workers()[p:length(workers())]
      @spawnat p2 begin
        A_on_p = @fetchfrom p A[:l]
        A[:l] = (I - A_on_p * A_on_p') * A[:l]
      end
    end
  end
end

# main
U = drand((m,m), workers(), [1,4])
orthonormalize(U)

```

See the line 4, instead I would like to perform the inplace orthonormalization of `A[:l]` (the local part of the input matrix).

If I use `@spawnat p LinearAlgebra.qr!(A[:l])` on line 4 instead, I don’t understand how to use `A[:l]` afterwards to orthogonalize the other partitions, because it contains only the reflectors. I cannot do `A_on_p * A_on_p'` for example.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [January 6, 2021, 10:16am UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/6 "2021-01-06T10:16:54Z")

</div>

You are correct that GS is slower than the Householder reflection QR in qr! GS is not the way to compute QR in general. However, recovering Q from the output of qr! is expensive.

Classical GS twice has very good parallel performance and is stable, more stable that the variants of modified GS in the codes. GMRES in Trilinos, for example, is done that way,

---

<div class="post-metadata">

**Author:** ![matthiasbe](https://avatars.discourse-cdn.com/v4/letter/m/5fc32e/32.png) [@matthiasbe](https://discourse.julialang.org/u/matthiasbe)\
**Post date:** [January 11, 2021, 3:57pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/7 "2021-01-11T15:57:59Z")

</div>

I found out that projecting when using the QR compressed format ([`QRCompactWY`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.QRCompactWY)) was never ending, rather than using the full matrix format was working

```julia
julia> using LinearAlgebra
julia> A = rand(1000,100)
julia> fact = qr!(A)
julia> (I - fact.Q * fact.Q') # Computing the projector in this way is very long (I didn't see the end)
julia> Q = Matrix(fact.Q)
julia> I - Q * Q' # In a matrix form, it takes a few ms

```

I would have expected using the compressed format to be faster, because it use householders transformations. I don’t see yet why it would be so long in practice.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 11, 2021, 4:34pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/8 "2021-01-11T16:34:11Z")

</div>

> [@matthiasbe](#):
>
> I would have expected using the compressed format to be faster, because it use householders transformations. I don’t see yet why it would be so long in practice.

The compressed format is designed for the circumstance where you don’t try to explicitly compute the matrix. For example, the projector `P = I - fact.Q * fact.Q'` should never be explicitly computed — rather, `P*x` should be computed implicitly as `x - fact.Q * (fact.Q'*x)` or similar.

---

<div class="post-metadata">

**Author:** ![matthiasbe](https://avatars.discourse-cdn.com/v4/letter/m/5fc32e/32.png) [@matthiasbe](https://discourse.julialang.org/u/matthiasbe)\
**Post date:** [January 11, 2021, 4:48pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/9 "2021-01-11T16:48:40Z")

</div>

@stevengj thank you for this remark, it makes sens indeed.

Using parenthesis like you pointed out solved the problem of undefinite running.

This led me to another problem, now I can compute the result of the projection with both formats, but the results do not concord

```julia
julia> using LinearAlgebra
julia> A = rand(1000,100)
julia> fact = qr!(A)
julia> B = rand(1000,100)
julia> B1 = B - fact.Q * (fact.Q' * B) # B1 uses the QRCompactWY format
julia> Q = Matrix(fact.Q)
julia> B2 = B - Q * (Q' * B) # B2 uses the explicit matrix format
julia> norm(B1 - B2, Inf)
0.8313667296907561 # The output shows both give very different results

```

I check the results in B1 and B2, and only B2 contains the expected result

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 11, 2021, 5:09pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/10 "2021-01-11T17:09:18Z")

</div>

> [@matthiasbe](#):
>
> I check the results in B1 and B2, and only B1 contains the expected result

I think you mean only B2 contains the expected result?

I forgot that `fact.Q' * B` by default computes the product of the “full” (square) Q^\* matrix with B — so that QQ^\* B = B — whereas you want the “thin” Q matrix result \hat{Q}^\* B, corresponding to the first 100 rows of Q^\* B:

```julia
B3 = B - fact.Q * (fact.Q' * B)[axes(B,2),:];

```

(I’m not sure offhand if there is a more efficient way to get this?)

---

<div class="post-metadata">

**Author:** ![matthiasbe](https://avatars.discourse-cdn.com/v4/letter/m/5fc32e/32.png) [@matthiasbe](https://discourse.julialang.org/u/matthiasbe)\
**Post date:** [January 11, 2021, 6:21pm UTC](https://discourse.julialang.org/t/orthonormalize-a-matrix-in-place/52887/11 "2021-01-11T18:21:32Z")

</div>

Exactly I meant B2, sorry.

Ok this explain the difference, thank you !

For information I had a speedup of circa 2 to 3 with using the compact format.