]>
gitweb.michael.orlitzky.com - sage.d.git/blob - eja_utils.py
7c6c581329dd64aaa8f3d958e00b2f9eb3221fcb
1 from sage
.functions
.other
import sqrt
2 from sage
.matrix
.constructor
import matrix
3 from sage
.modules
.free_module_element
import vector
7 Scale the vector, matrix, or cartesian-product-of-those-things
10 This works around the inability to scale certain elements of
11 Cartesian product spaces, as reported in
13 https://trac.sagemath.org/ticket/31435
17 This will do the wrong thing if you feed it a tuple or list.
21 sage: from mjo.eja.eja_utils import _scale
25 sage: v = vector(QQ, (1,2,3))
28 sage: m = matrix(QQ, [[1,2],[3,4]])
29 sage: M = cartesian_product([m.parent(), m.parent()])
30 sage: _scale(M((m,m)), 2)
36 if hasattr(x
, 'cartesian_factors'):
38 return P(tuple( _scale(x_i
, alpha
)
39 for x_i
in x
.cartesian_factors() ))
46 Flatten a vector, matrix, or cartesian product of those things
49 If the entries of the matrix themselves belong to a real vector
50 space (such as the complex numbers which can be thought of as
51 pairs of real numbers), they will also be expanded in vector form
52 and flattened into the list.
56 sage: from mjo.eja.eja_utils import _all2list
57 sage: from mjo.hurwitz import (QuaternionMatrixAlgebra,
59 ....: OctonionMatrixAlgebra)
63 sage: _all2list([[1]])
68 sage: V1 = VectorSpace(QQ,2)
69 sage: V2 = MatrixSpace(QQ,2)
73 sage: y2 = V2([0,1,1,0])
74 sage: _all2list((x1,y1))
76 sage: _all2list((x2,y2))
78 sage: M = cartesian_product([V1,V2])
79 sage: _all2list(M((x1,y1)))
81 sage: _all2list(M((x2,y2)))
86 sage: _all2list(Octonions().one())
87 [1, 0, 0, 0, 0, 0, 0, 0]
88 sage: _all2list(OctonionMatrixAlgebra(1).one())
89 [1, 0, 0, 0, 0, 0, 0, 0]
93 sage: _all2list(QuaternionAlgebra(QQ, -1, -1).one())
95 sage: _all2list(QuaternionMatrixAlgebra(1).one())
100 sage: V1 = VectorSpace(QQ,2)
101 sage: V2 = OctonionMatrixAlgebra(1,field=QQ)
102 sage: C = cartesian_product([V1,V2])
105 sage: _all2list(C( (x1,y1) ))
106 [3, 4, 1, 0, 0, 0, 0, 0, 0, 0]
109 if hasattr(x
, 'to_vector'):
110 # This works on matrices of e.g. octonions directly, without
111 # first needing to convert them to a list of octonions and
112 # then recursing down into the list. It also avoids the wonky
113 # list(x) when x is an element of a CFM. I don't know what it
114 # returns but it aint the coordinates. This will fall through
115 # to the iterable case the next time around.
116 return _all2list(x
.to_vector())
120 except TypeError: # x is not iterable
124 # Avoid the retardation of list(QQ(1)) == [1].
127 return sum(list( map(_all2list
, xl
) ), [])
132 return vector(m
.base_ring(), m
.list())
135 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
137 def gram_schmidt(v
, inner_product
=None):
139 Perform Gram-Schmidt on the list ``v`` which are assumed to be
140 vectors over the same base ring. Returns a list of orthonormalized
141 vectors over the same base ring, which means that your base ring
142 needs to contain the appropriate roots.
146 sage: from mjo.eja.eja_utils import gram_schmidt
150 If you start with an orthonormal set, you get it back. We can use
151 the rationals here because we don't need any square roots::
153 sage: v1 = vector(QQ, (1,0,0))
154 sage: v2 = vector(QQ, (0,1,0))
155 sage: v3 = vector(QQ, (0,0,1))
157 sage: gram_schmidt(v) == v
160 The usual inner-product and norm are default::
162 sage: v1 = vector(AA,(1,2,3))
163 sage: v2 = vector(AA,(1,-1,6))
164 sage: v3 = vector(AA,(2,1,-1))
166 sage: u = gram_schmidt(v)
167 sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
169 sage: bool(u[0].inner_product(u[1]) == 0)
171 sage: bool(u[0].inner_product(u[2]) == 0)
173 sage: bool(u[1].inner_product(u[2]) == 0)
177 But if you supply a custom inner product, the result is
178 orthonormal with respect to that (and not the usual inner
181 sage: v1 = vector(AA,(1,2,3))
182 sage: v2 = vector(AA,(1,-1,6))
183 sage: v3 = vector(AA,(2,1,-1))
185 sage: B = matrix(AA, [ [6, 4, 2],
188 sage: ip = lambda x,y: (B*x).inner_product(y)
189 sage: norm = lambda x: ip(x,x)
190 sage: u = gram_schmidt(v,ip)
191 sage: all( norm(u_i) == 1 for u_i in u )
193 sage: ip(u[0],u[1]).is_zero()
195 sage: ip(u[0],u[2]).is_zero()
197 sage: ip(u[1],u[2]).is_zero()
200 This Gram-Schmidt routine can be used on matrices as well, so long
201 as an appropriate inner-product is provided::
203 sage: E11 = matrix(AA, [ [1,0],
205 sage: E12 = matrix(AA, [ [0,1],
207 sage: E22 = matrix(AA, [ [0,0],
209 sage: I = matrix.identity(AA,2)
210 sage: trace_ip = lambda X,Y: (X*Y).trace()
211 sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
213 [1 0] [ 0 0.7071067811865475?] [0 0]
214 [0 0], [0.7071067811865475? 0], [0 1]
217 It even works on Cartesian product spaces whose factors are vector
220 sage: V1 = VectorSpace(AA,2)
221 sage: V2 = MatrixSpace(AA,2)
222 sage: M = cartesian_product([V1,V2])
224 sage: x2 = V1([1,-1])
226 sage: y2 = V2([0,1,1,0])
227 sage: z1 = M((x1,y1))
228 sage: z2 = M((x2,y2))
230 ....: return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
231 sage: U = gram_schmidt([z1,z2], inner_product=ip)
241 Ensure that zero vectors don't get in the way::
243 sage: v1 = vector(AA,(1,2,3))
244 sage: v2 = vector(AA,(1,-1,6))
245 sage: v3 = vector(AA,(0,0,0))
247 sage: len(gram_schmidt(v)) == 2
250 if inner_product
is None:
251 inner_product
= lambda x
,y
: x
.inner_product(y
)
253 ip
= inner_product(x
,x
)
254 # Don't expand the given field; the inner-product's codomain
255 # is already correct. For example QQ(2).sqrt() returns sqrt(2)
256 # in SR, and that will give you weird errors about symbolics
257 # when what's really going wrong is that you're trying to
258 # orthonormalize in QQ.
259 return ip
.parent()(ip
.sqrt())
261 v
= list(v
) # make a copy, don't clobber the input
263 # Drop all zero vectors before we start.
264 v
= [ v_i
for v_i
in v
if not v_i
.is_zero() ]
272 # Our "zero" needs to belong to the right space for sum() to work.
273 zero
= v
[0].parent().zero()
276 if hasattr(v
[0], 'cartesian_factors'):
277 # Only use the slow implementation if necessary.
281 return sc(x
, (inner_product(x
,y
)/inner_product(x
,x
)))
283 # First orthogonalize...
284 for i
in range(1,len(v
)):
285 # Earlier vectors can be made into zero so we have to ignore them.
286 v
[i
] -= sum( (proj(v
[j
],v
[i
])
288 if not v
[j
].is_zero() ),
291 # And now drop all zero vectors again if they were "orthogonalized out."
292 v
= [ v_i
for v_i
in v
if not v_i
.is_zero() ]
294 # Just normalize. If the algebra is missing the roots, we can't add
295 # them here because then our subalgebra would have a bigger field
296 # than the superalgebra.
297 for i
in range(len(v
)):
298 v
[i
] = sc(v
[i
], ~
norm(v
[i
]))