Gram Schmidit python to julia

General Usage Performance
1

Hello, good evening,

I’m having troubles to port an algorithm(Gram Schmidit for qr decompoosition) from python to julia. My results are going to infinity in Orthogonal Q matrix. May someone help?

	# create the matrix
	m = n = 4
	A = randn(m, n)

	Q = zeros(m, n)

	# the gs algorithm
	for i in 1:n
		for j in 1:i-1
			Q[:, i] = Q[:, i] .- dot(A[:, i], Q[:, j]) ./ dot(Q[:, j], Q[:, j]) .* Q[:, j]
		end
		# normalize
		Q[:, i] = Q[:, i] ./ norm(Q[:, i])
	end

	@info Qa
2

Posting the code as an image makes it quite a bit harder to replicate the results without typing everything in to test, so I might have missed something. But one problem that stands out is that you need to orthogonalize only against vectors that you have previously fully orthogonalized. I think changing your inner loop to

for j in 1:i-1

should get you closer and will likely fix the problem unless I did miss something. Proof reading code without running it is a challenge.

Edited because I didn’t proof read what I first wrote all that well…

1 Like
3

chaged it here but its going to nan yet. Some things in python seems to be very different from julia.

image
image722×594 41 KB

4

The edit to get rid of the image helped. The other problem that I missed is that you aren’t initializing Q correctly. I was focusing on the loops. Try

# create the matrix
m = n = 4
A = randn(m, n)

Q = A

# the gs algorithm
for i in 1:n
  for j in 1:i-1
    Q[:, i] = Q[:, i] .- dot(A[:, i], Q[:, j]) ./ dot(Q[:, j], Q[:, j]) .* Q[:, j]
  end
  # normalize
  Q[:, i] = Q[:, i] ./ norm(Q[:, i])
end

@info Q
1 Like
5

I was trying also here and found a similar solution also

image
image752×897 58.4 KB

6

But its a lot unstable. If you keep running the cell sometimes we get wrong results, but seems the issue is really initialization

image
image713×411 26.5 KB

7

I’m probably signing off for the night but happy to continue in the morning if it is helpful. I’d look at the signs of the columns of Q. A real QR factorization is only unique up to multiplication of the columns of Q by \pm 1. There is something to be said for deciding that the signs of the diagonal elements of R should be positive, which makes the entire decomposition unique for a matrix with linearly independent columns, but LAPACK does not guarantee this. And a basic Gram-Schmidt decomposition as you have here does. So there’s no reason to think Q .- Q_lib should be zero in every column. It will be zero in some columns where the signs on the diagonal of R happen to coincide and not in others. I think Python calls LAPACK, so there should be similar issues there using Python’s QR, although I’ve never really used Python seriously.

It’s also worth noting that Gram-Schmidt really is unstable and should only be used in specific circumstances, usually with reorthogonalization as needed to enforce numerical orthogonality. But that’s not going to be the issue here. That is a problem when A is ill-conditioned, which is very unlikely for the small random matrices you are generating.

1 Like
8

There are other implementations for this procedure in internet. I was porting python because its in a book I’m reading as its hard to find math content books in julia. Then I read python and write julia code.

9

I’d rather see Julia, but Python does seem now to be a common default over pseudocode for presenting algorithms. I’m curious: Was the loop error in the book? For me, I’m sometimes hesitant to tie the presentation of an algorithm to a specific language, but one of the advantages of using a real language instead of pseudocode is that you can run it and see if it works. However I’ve seen a lot of typos in code published in print. I have an impression that publishing code that was never run is fairly common.

1 Like
10

You are correct. Your code works because of the explanation you write. If you sum Q with Qstdlib, comparing the stdlib qr from julia, the nonzeros colums are zeros in another.

image
image713×415 28.5 KB

here is the python code


# create the matrix 
m = 4
n = 4
A = np.random.randn(m,n)

# initialize
Q = np.zeros((m,n))


# the GS algo
for i in range(n):
    
    # initialize
    Q[:,i] = A[:,i]
    
    # orthogonalize
    a = A[:,i] # convenience
    for j in range(i): # only to earlier cols
        q = Q[:,j] # convenience
        Q[:,i]=Q[:,i]-np.dot(a,q)/np.dot(q,q)*q
    
    # normalize
    Q[:,i] = Q[:,i] / np.linalg.norm(Q[:,i])
11

I see. I forgot about Python ranges. Those are fairly error prone to translate to Julia. I teach numerical analysis classes in which I look at and usually try to debug student code in Python, even if I don’t use it myself. I’m pretty good about switching between 1 based and 0 based indexing of arrays, but I still don’t spot errors with things like range(i) very easily.

1 Like