These computations are easy in Sage, with some quirks.
Dot products
v=vector([1,2,3])
w = vector([1,1,1])
v.dot_product(w)
v*w
v.norm()
Testing whether sets are orthogonal, orthonormal
The easiest way to test whether a set of vectors v_1, …, v_k is orthogonal is:
- Create a matrix A=[v_1 … v_k] with the vectors as its columns.
- Compute A^T*A, the product of the transpose of A with A. The result is diagonal if and only if v_1, …, v_k are orthogonal. The result is the identity if and only if v_1, …, v_k are orthonormal.
This works because the (i,j) entry of A^T*A is the dot product of v_i and v_j.
For example:
v1 = vector([1,1,1,1])
v2 = vector([1,-1,0,0])
v3 = vector([0,0,1,-1])
A = matrix([v1,v2,v3]).transpose()
A
A.transpose()*A
So, these vectors are orthogonal but not orthonormal (which is true–look at them).
v1 = vector([1/2,1/2,1/2,1/2])
v2 = vector([1/sqrt(2),-1/sqrt(2),0,0])
v3 = vector([0,0,1/sqrt(2),-1/sqrt(2)])
A = matrix([v1,v2,v3]).transpose()
A
A.transpose()*A
So, these vectors are orthonormal (again, obviously true).
v1 = vector([1,1,1,1])
v2 = vector([1,1,0,0])
v3 = vector([0,0,1,1])
A = matrix([v1,v2,v3]).transpose()
A
A.transpose()*A
So, these vectors are not orthogonal.
Projections
As far as I know, Sage does not have built in functions for finding orthogonal projects, but writing such functions is fairly easy. This is an optional homework problem.
Gram-Schmidt
Given a matrix A, “A.gram_schmidt(orthonormal=True)” applies the Gram-Schmidt process to the rows of A. The output is a pair of matrices (G,M), so that G is the result of the Gram-Schmidt process and A = M*G. We’re only interested in G in this course.
A second wrinkle is that the Gram-Schmidt process involves taking square roots. If you type a matrix with rational entries, Sage assumes you want to stay in the world of matrices with rational entries. So, for instance,
matrix([[1,2],[3,4]]).gram_schmidt(orthonormal=True)
gives an error.
There seem to be roughly two options to deal with this:
- Tell Sage you want to work with decimals. The result will be a matrix with ugly decimals, most likely.
- Tell Sage it’s okay if the resulting matrix has orthogonal but not orthonormal columns (or rather, rows). Then make the rows have length 1 on your own — it’s easy once they’re orthogonal.
For the first option, create your matrix with a command like
A = matrix(RDF,[[1,2],[3,4]])
(RDF stands for “real double field” — the collection of double-precision decimals.) Then
A.gram_schmidt(orthonormal=True)
will behave as expected.
For the second option, just leave off the “orthonormal=True”. For example:
A = matrix([[1,2],[3,4]])
A.gram_schmidt()
Again, note that Sage is applying Gram-Schmidt to the rows. The rows of the first matrix in the output are the result. Notice that they are orthogonal (but not length 1).
Finally, since that was all a bit of a pain, here’s some code which will do Gram-Schmidt to a list of vectors:
def better_gs(vectors):
A = matrix(vectors)
G = A.gram_schmidt()[0]
return [v / v.norm() for v in G.rows()]
(The indentation is important.) If you paste that into your notebook, applying better_gs to a list of vectors should have the result you expect:
v1 = vector([1,2])
v2 = vector([3,4])
better_gs([v1,v2])
