Julia Community 🟣

Vinod V
Vinod V

Posted on

Gram–Schmidt Orthogonalization with Julia

Theorem : Gram–Schmidt Let WW be a pp -dimensional subspace of RnR^n , and let w1,w2,…,wp{w_1, w_2, \dots , w_p} be any basis for WW . Then the set of vectors u1,u2,...,up{u_1, u_2, . . . , u_p} is an orthogonal basis for WW , where

u1=w1u_1 = w_1

u2=w2−u1Tw2u1Tu1u1u_2 = w_2 - \frac{u_1^Tw_2}{u_1^Tu_1}u_1

u3=w3−u1Tw3u1Tu1u1−u2Tw3u2Tu2u2u_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

ui=wi−∑k=1k=i−1ukTwiukTukuk,2≤i≤pu_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 wiw_i ,

wi=ui−∑k=1k=i−1ukTwiukTukuk,2≤i≤pw_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

[w1T⋮wnT]=[100⋯0t2110⋯0t31t321⋯0⋮⋮⋮⋮tn1tn2tn3⋯1][u1T⋮unT]\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}

tij=ujTwi−1/(ujTuj),∀(i>j);tij=1,∀(i=j);tij=0,∀(i<j)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,:]))
Enter fullscreen mode Exit fullscreen mode

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()
Enter fullscreen mode Exit fullscreen mode

Top comments (0)