Gram–Schmidt Orthogonalization with Julia

Theorem : Gram–Schmidt Let $W$ be a $p$ -dimensional subspace of $R^n$ , and let ${w_1, w_2, \dots , w_p}$ be any basis for $W$ . Then the set of vectors ${u_1, u_2, . . . , u_p}$ is an orthogonal basis for $W$ , where

$u_1 = w_1$

$u_2 = w_2 - \frac{u_1^Tw_2}{u_1^Tu_1}u_1$

$u_3 = w_3 - \frac{u_1^Tw_3}{u_1^Tu_1}u_1 - \frac{u_2^Tw_3}{u_2^Tu_2}u_2$

and where in general

$u_i = w_i - \sum_{k=1}^{k = i-1}\frac{u_k^Tw_i}{u_k^Tu_k}u_k, 2 \le i \le p$

By rewriting the above equations for $w_i$ ,

$w_i = u_i - \sum_{k=1}^{k = i-1}\frac{u_k^Tw_i}{u_k^Tu_k}u_k, 2 \le i \le p$

In matrix form

$\begin{bmatrix} w_1^T \\ \vdots \\ w_n^T\end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ t_{21} & 1 & 0 & \cdots & 0 \\ t_{31} & t_{32} & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots \\ t_{n1} & t_{n2} & t_{n3} & \cdots & 1 \end{bmatrix} \begin{bmatrix} u_1^T \\ \vdots \\ u_n^T\end{bmatrix}$

$t_{ij} = u_j^T w_{i-1}/(u_j^T u_j), \forall (i > j); t_{ij} = 1,\forall (i = j);t_{ij} = 0, \forall (i < j)$

Given below is an implementation of Gram-Schmidt orthogonalization in Julia

using LinearAlgebra

N = 10;
A = rand(N,N); #row vectors to be orthogonalized
U = zeros(N,N);
U[1,:] = A[1,:]; 
V = zeros(N,N); # for the for loop
V[1,:] = A[1,:]; 

U[2,:] = A[2,:] - dot(U[1,:],A[2,:])*U[1,:]/dot(U[1,:],U[1,:])

U[3,:] = A[3,:] - dot(U[1,:],A[3,:])*U[1,:]/dot(U[1,:],U[1,:]) - dot(U[2,:],A[3,:])*U[2,:]/dot(U[2,:],U[2,:])

println(dot(U[1,:],U[2,:]))

println(dot(U[1,:],U[3,:]))

println(dot(U[2,:],U[3,:]))

#for loop for a matrix of size N
for i in 2:N
    V[i,:] = A[i,:]
    for j in 1:(i-1)
        V[i,:] -= dot(V[j,:],A[i,:])*V[j,:]/dot(V[j,:],V[j,:])
    end
end


println(dot(V[1,:],V[2,:]))

println(dot(V[1,:],V[3,:]))

println(dot(V[2,:],V[3,:]))

println(dot(V[N,:],V[N-1,:]))

Using row operations

using LinearAlgebra

function test()

N = 4;A = rand(N,N);B = Matrix(1.0I,N,N);U = zeros(N,N);U[1,:] = A[1,:]; V = zeros(N,N);V[1,:] = A[1,:];

row_operations = [] 

B[2,1] = -dot(U[1,:],A[2,:])/dot(U[1,:],U[1,:])

push!(row_operations,B)

X = B*A

B = Matrix(1.0I,N,N)

B[3,1] = -dot(X[1,:],A[3,:])/dot(X[1,:],X[1,:])

B[3,2] = -dot(X[2,:],A[3,:])/dot(X[2,:],X[2,:])

push!(row_operations,B)

X = B*X

println(dot(X[1,:],X[2,:]))

println(dot(X[1,:],X[3,:]))

println(dot(X[2,:],X[3,:]))

X = zeros(N,N);X = copy(A);#X[1,:] = A[1,:];
row_operations = []
for i in 2:N
    C = Matrix(1.0I,N,N)
    for j in 1:(i-1)
        C[i,j] = -dot(X[j,:],A[i,:])/dot(X[j,:],X[j,:])     
    end
    push!(row_operations,C)
    X = C*X

end

println(dot(X[1,:],X[2,:]))

println(dot(X[1,:],X[3,:]))

println(dot(X[2,:],X[3,:]))

println(dot(X[N,:],X[N-1,:]))

end

test()

Julia Community 🟣

Gram–Schmidt Orthogonalization with Julia

Top comments (0)