implementing Gram-Schmidt process using SageMaths

179 Views Asked by At

I'm trying to implement a function myGramSchmidt(L), which takes a list L of vectors living in some inner product space, and returns a new list which has implemented the Gram-Schmidt process above.

(UPDATE) My code:

def myGramSchmidt(L):
    
    n = len(L)
    V = L.copy()
    for j in range(n):
        for i in range(j):
            V[j] -= ip(V[i],V[j]) * V[i]
        V[j]= V[j]/V[j].norm() #normalised vector
        print(V[j])
    print(ip(V[0],V[0]))
    print(ip(V[1],V[1]))
    print(ip(V[2],V[2]))
    return V

B=Matrix(QQ,3,3,[1,0,0, 0,2,0, 0,0,3])
V3 = VectorSpace(QQ, 3, inner_product_matrix = B)

v1 = V3(vector(random_matrix(QQ,1,3)))
v2 = V3(vector(random_matrix(QQ,1,3)))
v3 = V3(vector(random_matrix(QQ,1,3)))


L = [v1,v2,v3]
M1=Matrix(QQ,3,3,L)
pretty_print(M1)
myGramSchmidt(L)


I have generated matrix: $$ \left[ \begin{array}{ccc} \frac{1}{2}&0&-1 \\ 2&-1&-1 \\ 0&1&0 \end{array} \right] $$ with output:

(1/5*sqrt(5), 0, -2/5*sqrt(5))
(1/15*sqrt(6), -1/6*sqrt(6), 11/30*sqrt(6))
(2/15, 2/3, 11/15)

13/5
139/50
63/25

As you can see, the vectors are normalised, however the inner product of two same vectors does not equal to 1?

1

There are 1 best solutions below

4
On

The expected output should be: a list of vectors satisfying v.inner_product(v) == 1 and v.inner_product(w) == 0 when v and w are distinct.

You should interchange lines in your code: change

    V[j]= V[j]/V[j].norm() #normalised vector
    for i in range(j):
        V[j] -= ip(V[i],V[j]) * V[i] 

to

    for i in range(j):
        V[j] -= ip(V[i],V[j]) * V[i]
     V[j]= V[j]/V[j].norm() #normalised vector

(Don't normalize the vector until you've done all other modifications to it.)

Edit: a further problem is the norm method, which seems to ignore the inner product matrix. Try this instead:

    for i in range(j):
        V[j] -= ip(V[i],V[j]) * V[i]
     V[j]= V[j]/sqrt(V[j].inner_product(V[j])) #normalised vector

You can also add this as a check:

M = matrix(n, [V[i].inner_product(V[j]) for i in range(n) for j in range(n)])
print(M)