from sage.functions.other import sqrt from sage.matrix.constructor import matrix from sage.modules.free_module_element import vector def _scale(x, alpha): r""" Scale the vector, matrix, or cartesian-product-of-those-things ``x`` by ``alpha``. This works around the inability to scale certain elements of Cartesian product spaces, as reported in https://trac.sagemath.org/ticket/31435 ..WARNING: This will do the wrong thing if you feed it a tuple or list. SETUP:: sage: from mjo.eja.eja_utils import _scale EXAMPLES:: sage: v = vector(QQ, (1,2,3)) sage: _scale(v,2) (2, 4, 6) sage: m = matrix(QQ, [[1,2],[3,4]]) sage: M = cartesian_product([m.parent(), m.parent()]) sage: _scale(M((m,m)), 2) ([2 4] [6 8], [2 4] [6 8]) """ if hasattr(x, 'cartesian_factors'): P = x.parent() return P(tuple( _scale(x_i, alpha) for x_i in x.cartesian_factors() )) else: return x*alpha def _all2list(x): r""" Flatten a vector, matrix, or cartesian product of those things into a long list. If the entries of the matrix themselves belong to a real vector space (such as the complex numbers which can be thought of as pairs of real numbers), they will also be expanded in vector form and flattened into the list. SETUP:: sage: from mjo.eja.eja_utils import _all2list sage: from mjo.hurwitz import (QuaternionMatrixAlgebra, ....: Octonions, ....: OctonionMatrixAlgebra) EXAMPLES:: sage: _all2list([[1]]) [1] :: sage: V1 = VectorSpace(QQ,2) sage: V2 = MatrixSpace(QQ,2) sage: x1 = V1([1,1]) sage: x2 = V1([1,-1]) sage: y1 = V2.one() sage: y2 = V2([0,1,1,0]) sage: _all2list((x1,y1)) [1, 1, 1, 0, 0, 1] sage: _all2list((x2,y2)) [1, -1, 0, 1, 1, 0] sage: M = cartesian_product([V1,V2]) sage: _all2list(M((x1,y1))) [1, 1, 1, 0, 0, 1] sage: _all2list(M((x2,y2))) [1, -1, 0, 1, 1, 0] :: sage: _all2list(Octonions().one()) [1, 0, 0, 0, 0, 0, 0, 0] sage: _all2list(OctonionMatrixAlgebra(1).one()) [1, 0, 0, 0, 0, 0, 0, 0] :: sage: _all2list(QuaternionAlgebra(QQ, -1, -1).one()) [1, 0, 0, 0] sage: _all2list(QuaternionMatrixAlgebra(1).one()) [1, 0, 0, 0] :: sage: V1 = VectorSpace(QQ,2) sage: V2 = OctonionMatrixAlgebra(1,field=QQ) sage: C = cartesian_product([V1,V2]) sage: x1 = V1([3,4]) sage: y1 = V2.one() sage: _all2list(C( (x1,y1) )) [3, 4, 1, 0, 0, 0, 0, 0, 0, 0] """ if hasattr(x, 'to_vector'): # This works on matrices of e.g. octonions directly, without # first needing to convert them to a list of octonions and # then recursing down into the list. It also avoids the wonky # list(x) when x is an element of a CFM. I don't know what it # returns but it aint the coordinates. This will fall through # to the iterable case the next time around. return _all2list(x.to_vector()) try: xl = list(x) except TypeError: # x is not iterable return [x] if xl == [x]: # Avoid the retardation of list(QQ(1)) == [1]. return [x] return sum(list( map(_all2list, xl) ), []) def _mat2vec(m): return vector(m.base_ring(), m.list()) def _vec2mat(v): return matrix(v.base_ring(), sqrt(v.degree()), v.list()) 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 vectors over the same base ring, which means that your base ring needs to contain the appropriate roots. SETUP:: sage: from mjo.eja.eja_utils import gram_schmidt EXAMPLES: If you start with an orthonormal set, you get it back. We can use the rationals here because we don't need any square roots:: sage: v1 = vector(QQ, (1,0,0)) sage: v2 = vector(QQ, (0,1,0)) sage: v3 = vector(QQ, (0,0,1)) sage: v = [v1,v2,v3] sage: gram_schmidt(v) == v True The usual inner-product and norm are default:: sage: v1 = vector(AA,(1,2,3)) sage: v2 = vector(AA,(1,-1,6)) sage: v3 = vector(AA,(2,1,-1)) sage: v = [v1,v2,v3] sage: u = gram_schmidt(v) sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u ) True sage: bool(u[0].inner_product(u[1]) == 0) True sage: bool(u[0].inner_product(u[2]) == 0) True 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(AA,(1,2,3)) sage: v2 = vector(AA,(1,-1,6)) sage: v3 = vector(AA,(2,1,-1)) sage: v = [v1,v2,v3] sage: B = matrix(AA, [ [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 This Gram-Schmidt routine can be used on matrices as well, so long as an appropriate inner-product is provided:: sage: E11 = matrix(AA, [ [1,0], ....: [0,0] ]) sage: E12 = matrix(AA, [ [0,1], ....: [1,0] ]) sage: E22 = matrix(AA, [ [0,0], ....: [0,1] ]) sage: I = matrix.identity(AA,2) sage: trace_ip = lambda X,Y: (X*Y).trace() sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip) [ [1 0] [ 0 0.7071067811865475?] [0 0] [0 0], [0.7071067811865475? 0], [0 1] ] It even works on Cartesian product spaces whose factors are vector or matrix spaces:: sage: V1 = VectorSpace(AA,2) sage: V2 = MatrixSpace(AA,2) sage: M = cartesian_product([V1,V2]) sage: x1 = V1([1,1]) sage: x2 = V1([1,-1]) sage: y1 = V2.one() sage: y2 = V2([0,1,1,0]) sage: z1 = M((x1,y1)) sage: z2 = M((x2,y2)) sage: def ip(a,b): ....: return a[0].inner_product(b[0]) + (a[1]*b[1]).trace() sage: U = gram_schmidt([z1,z2], inner_product=ip) sage: ip(U[0],U[1]) 0 sage: ip(U[0],U[0]) 1 sage: ip(U[1],U[1]) 1 TESTS: Ensure that zero vectors don't get in the way:: sage: v1 = vector(AA,(1,2,3)) sage: v2 = vector(AA,(1,-1,6)) sage: v3 = vector(AA,(0,0,0)) sage: v = [v1,v2,v3] sage: len(gram_schmidt(v)) == 2 True """ if inner_product is None: inner_product = lambda x,y: x.inner_product(y) def norm(x): ip = inner_product(x,x) # Don't expand the given field; the inner-product's codomain # is already correct. For example QQ(2).sqrt() returns sqrt(2) # in SR, and that will give you weird errors about symbolics # when what's really going wrong is that you're trying to # orthonormalize in QQ. return ip.parent()(ip.sqrt()) v = list(v) # make a copy, don't clobber the input # Drop all zero vectors before we start. v = [ v_i for v_i in v if not v_i.is_zero() ] if len(v) == 0: # cool return v R = v[0].base_ring() # Our "zero" needs to belong to the right space for sum() to work. zero = v[0].parent().zero() sc = lambda x,a: a*x if hasattr(v[0], 'cartesian_factors'): # Only use the slow implementation if necessary. sc = _scale def proj(x,y): return sc(x, (inner_product(x,y)/inner_product(x,x))) # First orthogonalize... 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() ), zero ) # And now drop all zero vectors again if they were "orthogonalized out." v = [ v_i for v_i in v if not v_i.is_zero() ] # 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 range(len(v)): v[i] = sc(v[i], ~norm(v[i])) return v