]>
gitweb.michael.orlitzky.com - sage.d.git/blob - eja/eja_utils.py
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.octonions import Octonions, OctonionMatrixAlgebra
61 sage: _all2list([[1]])
66 sage: V1 = VectorSpace(QQ,2)
67 sage: V2 = MatrixSpace(QQ,2)
71 sage: y2 = V2([0,1,1,0])
72 sage: _all2list((x1,y1))
74 sage: _all2list((x2,y2))
76 sage: M = cartesian_product([V1,V2])
77 sage: _all2list(M((x1,y1)))
79 sage: _all2list(M((x2,y2)))
84 sage: _all2list(Octonions().one())
85 [1, 0, 0, 0, 0, 0, 0, 0]
86 sage: _all2list(OctonionMatrixAlgebra(1).one())
87 [1, 0, 0, 0, 0, 0, 0, 0]
91 sage: V1 = VectorSpace(QQ,2)
92 sage: V2 = OctonionMatrixAlgebra(1,field=QQ)
93 sage: C = cartesian_product([V1,V2])
96 sage: _all2list(C( (x1,y1) ))
97 [3, 4, 1, 0, 0, 0, 0, 0, 0, 0]
100 if hasattr(x
, 'to_vector'):
101 # This works on matrices of e.g. octonions directly, without
102 # first needing to convert them to a list of octonions and
103 # then recursing down into the list. It also avoids the wonky
104 # list(x) when x is an element of a CFM. I don't know what it
105 # returns but it aint the coordinates. This will fall through
106 # to the iterable case the next time around.
107 return _all2list(x
.to_vector())
111 except TypeError: # x is not iterable
115 # Avoid the retardation of list(QQ(1)) == [1].
118 return sum(list( map(_all2list
, xl
) ), [])
123 return vector(m
.base_ring(), m
.list())
126 return matrix(v
.base_ring(), sqrt(v
.degree()), v
.list())
128 def gram_schmidt(v
, inner_product
=None):
130 Perform Gram-Schmidt on the list ``v`` which are assumed to be
131 vectors over the same base ring. Returns a list of orthonormalized
132 vectors over the same base ring, which means that your base ring
133 needs to contain the appropriate roots.
137 sage: from mjo.eja.eja_utils import gram_schmidt
141 If you start with an orthonormal set, you get it back. We can use
142 the rationals here because we don't need any square roots::
144 sage: v1 = vector(QQ, (1,0,0))
145 sage: v2 = vector(QQ, (0,1,0))
146 sage: v3 = vector(QQ, (0,0,1))
148 sage: gram_schmidt(v) == v
151 The usual inner-product and norm are default::
153 sage: v1 = vector(AA,(1,2,3))
154 sage: v2 = vector(AA,(1,-1,6))
155 sage: v3 = vector(AA,(2,1,-1))
157 sage: u = gram_schmidt(v)
158 sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
160 sage: bool(u[0].inner_product(u[1]) == 0)
162 sage: bool(u[0].inner_product(u[2]) == 0)
164 sage: bool(u[1].inner_product(u[2]) == 0)
168 But if you supply a custom inner product, the result is
169 orthonormal with respect to that (and not the usual inner
172 sage: v1 = vector(AA,(1,2,3))
173 sage: v2 = vector(AA,(1,-1,6))
174 sage: v3 = vector(AA,(2,1,-1))
176 sage: B = matrix(AA, [ [6, 4, 2],
179 sage: ip = lambda x,y: (B*x).inner_product(y)
180 sage: norm = lambda x: ip(x,x)
181 sage: u = gram_schmidt(v,ip)
182 sage: all( norm(u_i) == 1 for u_i in u )
184 sage: ip(u[0],u[1]).is_zero()
186 sage: ip(u[0],u[2]).is_zero()
188 sage: ip(u[1],u[2]).is_zero()
191 This Gram-Schmidt routine can be used on matrices as well, so long
192 as an appropriate inner-product is provided::
194 sage: E11 = matrix(AA, [ [1,0],
196 sage: E12 = matrix(AA, [ [0,1],
198 sage: E22 = matrix(AA, [ [0,0],
200 sage: I = matrix.identity(AA,2)
201 sage: trace_ip = lambda X,Y: (X*Y).trace()
202 sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
204 [1 0] [ 0 0.7071067811865475?] [0 0]
205 [0 0], [0.7071067811865475? 0], [0 1]
208 It even works on Cartesian product spaces whose factors are vector
211 sage: V1 = VectorSpace(AA,2)
212 sage: V2 = MatrixSpace(AA,2)
213 sage: M = cartesian_product([V1,V2])
215 sage: x2 = V1([1,-1])
217 sage: y2 = V2([0,1,1,0])
218 sage: z1 = M((x1,y1))
219 sage: z2 = M((x2,y2))
221 ....: return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
222 sage: U = gram_schmidt([z1,z2], inner_product=ip)
232 Ensure that zero vectors don't get in the way::
234 sage: v1 = vector(AA,(1,2,3))
235 sage: v2 = vector(AA,(1,-1,6))
236 sage: v3 = vector(AA,(0,0,0))
238 sage: len(gram_schmidt(v)) == 2
241 if inner_product
is None:
242 inner_product
= lambda x
,y
: x
.inner_product(y
)
244 ip
= inner_product(x
,x
)
245 # Don't expand the given field; the inner-product's codomain
246 # is already correct. For example QQ(2).sqrt() returns sqrt(2)
247 # in SR, and that will give you weird errors about symbolics
248 # when what's really going wrong is that you're trying to
249 # orthonormalize in QQ.
250 return ip
.parent()(ip
.sqrt())
252 v
= list(v
) # make a copy, don't clobber the input
254 # Drop all zero vectors before we start.
255 v
= [ v_i
for v_i
in v
if not v_i
.is_zero() ]
263 # Our "zero" needs to belong to the right space for sum() to work.
264 zero
= v
[0].parent().zero()
267 if hasattr(v
[0], 'cartesian_factors'):
268 # Only use the slow implementation if necessary.
272 return sc(x
, (inner_product(x
,y
)/inner_product(x
,x
)))
274 # First orthogonalize...
275 for i
in range(1,len(v
)):
276 # Earlier vectors can be made into zero so we have to ignore them.
277 v
[i
] -= sum( (proj(v
[j
],v
[i
])
279 if not v
[j
].is_zero() ),
282 # And now drop all zero vectors again if they were "orthogonalized out."
283 v
= [ v_i
for v_i
in v
if not v_i
.is_zero() ]
285 # Just normalize. If the algebra is missing the roots, we can't add
286 # them here because then our subalgebra would have a bigger field
287 # than the superalgebra.
288 for i
in range(len(v
)):
289 v
[i
] = sc(v
[i
], ~
norm(v
[i
]))