]>
gitweb.michael.orlitzky.com - sage.d.git/blob - eja/eja_utils.py
1 from sage
.structure
.element
import is_Matrix
5 Scale the vector, matrix, or cartesian-product-of-those-things
8 This works around the inability to scale certain elements of
9 Cartesian product spaces, as reported in
11 https://trac.sagemath.org/ticket/31435
15 This will do the wrong thing if you feed it a tuple or list.
19 sage: from mjo.eja.eja_utils import _scale
23 sage: v = vector(QQ, (1,2,3))
26 sage: m = matrix(QQ, [[1,2],[3,4]])
27 sage: M = cartesian_product([m.parent(), m.parent()])
28 sage: _scale(M((m,m)), 2)
34 if hasattr(x
, 'cartesian_factors'):
36 return P(tuple( _scale(x_i
, alpha
)
37 for x_i
in x
.cartesian_factors() ))
44 Flatten a vector, matrix, or cartesian product of those things
47 If the entries of the matrix themselves belong to a real vector
48 space (such as the complex numbers which can be thought of as
49 pairs of real numbers), they will also be expanded in vector form
50 and flattened into the list.
54 sage: from mjo.eja.eja_utils import _all2list
55 sage: from mjo.hurwitz import (QuaternionMatrixAlgebra,
57 ....: 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: _all2list(QuaternionAlgebra(QQ, -1, -1).one())
93 sage: _all2list(QuaternionMatrixAlgebra(1).one())
98 sage: V1 = VectorSpace(QQ,2)
99 sage: V2 = OctonionMatrixAlgebra(1,field=QQ)
100 sage: C = cartesian_product([V1,V2])
103 sage: _all2list(C( (x1,y1) ))
104 [3, 4, 1, 0, 0, 0, 0, 0, 0, 0]
107 if hasattr(x
, 'to_vector'):
108 # This works on matrices of e.g. octonions directly, without
109 # first needing to convert them to a list of octonions and
110 # then recursing down into the list. It also avoids the wonky
111 # list(x) when x is an element of a CFM. I don't know what it
112 # returns but it aint the coordinates. We don't recurse
113 # because vectors can only contain ring elements as entries.
114 return x
.to_vector().list()
117 # This sucks, but for performance reasons we don't want to
118 # call _all2list recursively on the contents of a matrix
119 # when we don't have to (they only contain ring elements
125 except TypeError: # x is not iterable
129 # Avoid the retardation of list(QQ(1)) == [1].
132 return sum( map(_all2list
, xl
) , [])
135 def gram_schmidt(v
, inner_product
=None):
137 Perform Gram-Schmidt on the list ``v`` which are assumed to be
138 vectors over the same base ring. Returns a list of orthonormalized
139 vectors over the same base ring, which means that your base ring
140 needs to contain the appropriate roots.
144 sage: from mjo.eja.eja_utils import gram_schmidt
148 If you start with an orthonormal set, you get it back. We can use
149 the rationals here because we don't need any square roots::
151 sage: v1 = vector(QQ, (1,0,0))
152 sage: v2 = vector(QQ, (0,1,0))
153 sage: v3 = vector(QQ, (0,0,1))
155 sage: gram_schmidt(v) == v
158 The usual inner-product and norm are default::
160 sage: v1 = vector(AA,(1,2,3))
161 sage: v2 = vector(AA,(1,-1,6))
162 sage: v3 = vector(AA,(2,1,-1))
164 sage: u = gram_schmidt(v)
165 sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
167 sage: bool(u[0].inner_product(u[1]) == 0)
169 sage: bool(u[0].inner_product(u[2]) == 0)
171 sage: bool(u[1].inner_product(u[2]) == 0)
175 But if you supply a custom inner product, the result is
176 orthonormal with respect to that (and not the usual inner
179 sage: v1 = vector(AA,(1,2,3))
180 sage: v2 = vector(AA,(1,-1,6))
181 sage: v3 = vector(AA,(2,1,-1))
183 sage: B = matrix(AA, [ [6, 4, 2],
186 sage: ip = lambda x,y: (B*x).inner_product(y)
187 sage: norm = lambda x: ip(x,x)
188 sage: u = gram_schmidt(v,ip)
189 sage: all( norm(u_i) == 1 for u_i in u )
191 sage: ip(u[0],u[1]).is_zero()
193 sage: ip(u[0],u[2]).is_zero()
195 sage: ip(u[1],u[2]).is_zero()
198 This Gram-Schmidt routine can be used on matrices as well, so long
199 as an appropriate inner-product is provided::
201 sage: E11 = matrix(AA, [ [1,0],
203 sage: E12 = matrix(AA, [ [0,1],
205 sage: E22 = matrix(AA, [ [0,0],
207 sage: I = matrix.identity(AA,2)
208 sage: trace_ip = lambda X,Y: (X*Y).trace()
209 sage: gram_schmidt([E11,E12,I,E22], inner_product=trace_ip)
211 [1 0] [ 0 0.7071067811865475?] [0 0]
212 [0 0], [0.7071067811865475? 0], [0 1]
215 It even works on Cartesian product spaces whose factors are vector
218 sage: V1 = VectorSpace(AA,2)
219 sage: V2 = MatrixSpace(AA,2)
220 sage: M = cartesian_product([V1,V2])
222 sage: x2 = V1([1,-1])
224 sage: y2 = V2([0,1,1,0])
225 sage: z1 = M((x1,y1))
226 sage: z2 = M((x2,y2))
228 ....: return a[0].inner_product(b[0]) + (a[1]*b[1]).trace()
229 sage: U = gram_schmidt([z1,z2], inner_product=ip)
239 Ensure that zero vectors don't get in the way::
241 sage: v1 = vector(AA,(1,2,3))
242 sage: v2 = vector(AA,(1,-1,6))
243 sage: v3 = vector(AA,(0,0,0))
245 sage: len(gram_schmidt(v)) == 2
255 if inner_product
is None:
256 inner_product
= lambda x
,y
: x
.inner_product(y
)
259 if hasattr(V
, 'cartesian_factors'):
260 # Only use the slow implementation if necessary.
264 # project y onto the span of {x}
265 return sc(x
, (inner_product(x
,y
)/inner_product(x
,x
)))
268 # Don't extend the given field with the necessary
269 # square roots. This will probably throw weird
270 # errors about the symbolic ring if you e.g. try
271 # to use it on a set of rational vectors that isn't
272 # already orthonormalized.
273 return sc(x
, ~
inner_product(x
,x
).sqrt())
275 v_out
= [] # make a copy, don't clobber the input
277 for (i
, v_i
) in enumerate(v
):
278 ortho_v_i
= v_i
- V
.sum( proj(v_out
[j
],v_i
) for j
in range(i
) )
279 if not ortho_v_i
.is_zero():
280 v_out
.append(normalize(ortho_v_i
))