From 3c1c9170143a6412b050602cd79a383e3da6c821 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 27 Nov 2020 09:16:11 -0500 Subject: [PATCH] eja: allow non-standard inner product in gram_schmidt. --- mjo/eja/eja_utils.py | 42 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 36 insertions(+), 6 deletions(-) diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py index 49e3078..4d70e06 100644 --- a/mjo/eja/eja_utils.py +++ b/mjo/eja/eja_utils.py @@ -11,7 +11,7 @@ def _mat2vec(m): def _vec2mat(v): return matrix(v.base_ring(), sqrt(v.degree()), v.list()) -def gram_schmidt(v): +def gram_schmidt(v, inner_product=None): """ Perform Gram-Schmidt on the list ``v`` which are assumed to be vectors over the same base ring. Returns a list of orthonormalized @@ -22,7 +22,9 @@ def gram_schmidt(v): sage: from mjo.eja.eja_utils import gram_schmidt - EXAMPLES:: + EXAMPLES: + + The usual inner-product and norm are default:: sage: v1 = vector(QQ,(1,2,3)) sage: v2 = vector(QQ,(1,-1,6)) @@ -38,6 +40,30 @@ def gram_schmidt(v): sage: bool(u[1].inner_product(u[2]) == 0) True + + But if you supply a custom inner product, the result is + orthonormal with respect to that (and not the usual inner + product):: + + sage: v1 = vector(QQ,(1,2,3)) + sage: v2 = vector(QQ,(1,-1,6)) + sage: v3 = vector(QQ,(2,1,-1)) + sage: v = [v1,v2,v3] + sage: B = matrix(QQ, [ [6, 4, 2], + ....: [4, 5, 4], + ....: [2, 4, 9] ]) + sage: ip = lambda x,y: (B*x).inner_product(y) + sage: norm = lambda x: ip(x,x) + sage: u = gram_schmidt(v,ip) + sage: all( norm(u_i) == 1 for u_i in u ) + True + sage: ip(u[0],u[1]).is_zero() + True + sage: ip(u[0],u[2]).is_zero() + True + sage: ip(u[1],u[2]).is_zero() + True + TESTS: Ensure that zero vectors don't get in the way:: @@ -50,8 +76,12 @@ def gram_schmidt(v): True """ + if inner_product is None: + inner_product = lambda x,y: x.inner_product(y) + norm = lambda x: inner_product(x,x).sqrt() + def proj(x,y): - return (y.inner_product(x)/x.inner_product(x))*x + return (inner_product(x,y)/inner_product(x,x))*x v = list(v) # make a copy, don't clobber the input @@ -65,7 +95,7 @@ def gram_schmidt(v): R = v[0].base_ring() # First orthogonalize... - for i in xrange(1,len(v)): + for i in range(1,len(v)): # Earlier vectors can be made into zero so we have to ignore them. v[i] -= sum( proj(v[j],v[i]) for j in range(i) if not v[j].is_zero() ) @@ -75,7 +105,7 @@ def gram_schmidt(v): # Just normalize. If the algebra is missing the roots, we can't add # them here because then our subalgebra would have a bigger field # than the superalgebra. - for i in xrange(len(v)): - v[i] = v[i] / v[i].norm() + for i in range(len(v)): + v[i] = v[i] / norm(v[i]) return v -- 2.44.2