Eigen spaces of algebraic matrices broken ?
129 views
Skip to first unread message
Emmanuel Charpentier
Jan 29, 2022, 9:34:05 AM1/29/22
to sage-support
Setup : Sage 9.5.rc1 running in Debian testing on core i7 + 16 GB RAM.
def test(Size=2, Ring=QQ):
from time import time as stime
with seed(0):
M = matrix(Ring, Size, Size, lambda u, v:Ring.random_element())
t0 = stime()
SL = M.eigenspaces_left(algebraic_multiplicity=True)
t1 = stime()
VL = M.eigenvectors_left()
t2 = stime()
return t1-t0, SL, t2-t1, VL
test() runs as expected, as well as test(Size=5) and test(Ring=AA), with very reasonable runtimes (under a second, IIRC). But test(Size=3, Ring=AA) “never returns” (meaning that it hadn’t returned when I interrupted it after tens of minutes). Further exploration showec that neither M.eigenspaces_left(algebraic_multiplicity=True) nor M.eigenvectors_left() return.
Is that an expected behaviour ?
FWIW, stack trace at interruption :
sage: AA3=test(Size=3, Ring=AA)
AA3=test(Size=3, Ring=AA) C-c C-c---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in minpoly(self)
4492 try:
-> 4493 return self._minimal_polynomial
4494 except AttributeError:
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__getattr__ (build/cythonized/sage/structure/element.c:4754)()
493 """
--> 494 return self.getattr_from_category(name)
495
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.getattr_from_category (build/cythonized/sage/structure/element.c:4866)()
506 cls = P._abstract_element_class
--> 507 return getattr_from_other_class(self, cls, name)
508
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/cpython/getattr.pyx in sage.cpython.getattr.getattr_from_other_class (build/cythonized/sage/cpython/getattr.c:2633)()
360 dummy_error_message.name = name
--> 361 raise AttributeError(dummy_error_message)
362 attribute = <object>attr
AttributeError: 'sage.rings.real_mpfi.RealIntervalFieldElement' object has no attribute '_evaluate_polynomial'
During handling of the above exception, another exception occurred:
KeyboardInterrupt Traceback (most recent call last)
<ipython-input-344-ff7c46010ad2> in <module>
----> 1 AA3=test(Size=Integer(3), Ring=AA)
<string> in test(Size, Ring)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.eigenspaces_left (build/cythonized/sage/matrix/matrix2.c:43588)()
6270 self = self.change_ring(F)
6271 A = self - alpha
-> 6272 W = A.kernel()
6273 V.append((alpha, W.ambient_module().span_of_basis(W.basis()), e))
6274 else:
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.left_kernel (build/cythonized/sage/matrix/matrix2.c:32440)()
4947
4948 tm = verbose("computing left kernel for %sx%s matrix" % (self.nrows(), self.ncols()),level=1)
-> 4949 K = self.transpose().right_kernel(*args, **kwds)
4950 self.cache('left_kernel', K)
4951 verbose("done computing left kernel for %sx%s matrix" % (self.nrows(), self.ncols()),level=1,t=tm)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.right_kernel (build/cythonized/sage/matrix/matrix2.c:31909)()
4785
4786 # Go get the kernel matrix, this is where it all happens
-> 4787 M = self.right_kernel_matrix(*args, **kwds)
4788
4789 ambient = R**self.ncols()
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.right_kernel_matrix (build/cythonized/sage/matrix/matrix2.c:30829)()
4394
4395 if M is None and R in _Fields:
-> 4396 format, M = self._right_kernel_matrix_over_field()
4397
4398 if M is None and R.is_integral_domain():
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix._right_kernel_matrix_over_field (build/cythonized/sage/matrix/matrix2.c:27969)()
3686 from sage.matrix.matrix_space import MatrixSpace
3687 tm = verbose("computing right kernel matrix over an arbitrary field for %sx%s matrix" % (self.nrows(), self.ncols()),level=1)
-> 3688 E = self.echelon_form(*args, **kwds)
3689 pivots = E.pivots()
3690 pivots_set = set(pivots)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.echelon_form (build/cythonized/sage/matrix/matrix2.c:52440)()
7686 E = self.__copy__()
7687 if algorithm == 'default':
-> 7688 v = E.echelonize(cutoff=cutoff, **kwds)
7689 else:
7690 v = E.echelonize(algorithm = algorithm, cutoff=cutoff, **kwds)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.echelonize (build/cythonized/sage/matrix/matrix2.c:51817)()
7583 if self.base_ring() in _Fields:
7584 if algorithm in ['classical', 'partial_pivoting', 'scaled_partial_pivoting', 'scaled_partial_pivoting_valuation']:
-> 7585 self._echelon_in_place(algorithm)
7586 elif algorithm == 'strassen':
7587 self._echelon_strassen(cutoff)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix._echelon_in_place (build/cythonized/sage/matrix/matrix2.c:53308)()
7848 sig_check()
7849 for r in range(start_row, nr):
-> 7850 if A.get_unsafe(r, c):
7851 pivots.append(c)
7852 a_inverse = ~A.get_unsafe(r, c)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in __bool__(self)
4043 right = sd._right if op is operator.sub else -sd._right
4044
-> 4045 lp = left.minpoly()
4046 rp = right.minpoly()
4047 if lp != rp:
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in minpoly(self)
4493 return self._minimal_polynomial
4494 except AttributeError:
-> 4495 self.exactify()
4496 self._minimal_polynomial = self._descr.minpoly()
4497 return self._minimal_polynomial
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8430 right = self._right
8431 left.exactify()
-> 8432 right.exactify()
8433 gen = left._exact_field().union(right._exact_field())
8434 left_value = gen(left._exact_value())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8430 right = self._right
8431 left.exactify()
-> 8432 right.exactify()
8433 gen = left._exact_field().union(right._exact_field())
8434 left_value = gen(left._exact_value())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8430 right = self._right
8431 left.exactify()
-> 8432 right.exactify()
8433 gen = left._exact_field().union(right._exact_field())
8434 left_value = gen(left._exact_value())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8176
8177 if op == '~':
-> 8178 arg.exactify()
8179 return arg._descr.invert(None)
8180
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8431 left.exactify()
8432 right.exactify()
-> 8433 gen = left._exact_field().union(right._exact_field())
8434 left_value = gen(left._exact_value())
8435 right_value = gen(right._exact_value())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in union(self, other)
3230 op = QQx(op)
3231 # pari_nf = self._field.pari_nf()
-> 3232 pari_nf = self.pari_field()
3233 factors = list(pari_nf.nffactor(op).lift())[0]
3234 x, y = QQxy.gens()
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in pari_field(self)
3134 if self._pari_field is None:
3135 pari_pol = self._field.pari_polynomial("y")
-> 3136 self._pari_field = pari_pol.nfinit(1)
3137 return self._pari_field
3138
cypari2/auto_gen.pxi in cypari2.gen.Gen_base.nfinit()
KeyboardInterrupt:
Any hint appreciated…
Emmanuel Charpentier
Jan 29, 2022, 4:51:14 PM1/29/22
to sage-support
The promlem seems tolie wit (my use of) polynomial rings to compute the eigenvectors. Manually, after executing
dims = M.dimensions()
if dims[0] != dims[1]: raise DomainError("Not a square matrix !")
dim = dims[0]
BR = M.base_ring()
try:
WR = BR.algebraic_closure()
except AttributeError:
WR = BR
WPR = WR["t"]
t = WPR.gens()[0]
# CP = WPR(M.charpoly(var="t"))
# Eigenvalues
CP = M.charpoly(var="t")
lambdas = CP.roots(ring=WR)
# Eigenvectors
Id = identity_matrix(WR,dim)
WM = M.apply_map(lambda u:WR(u))
PR = PolynomialRing(WR,list(("v%d"%u for u in range(dim))))
J=PR.ideal(*(WM-Id*lambdas[0][0])*vector(PR.gens()))
neither J.dimension() nor J.variety() returns ; the interruption gives :
sage: J.dimension()
C-c C-c---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
<ipython-input-684-2e62bddf16b7> in <module>
----> 1 J.dimension()
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/polynomial/multi_polynomial_ideal.py in __call__(self, *args, **kwds)
295 if not R.base_ring().is_field():
296 raise ValueError("Coefficient ring must be a field for function '%s'."%(self.f.__name__))
--> 297 return self.f(self._instance, *args, **kwds)
298
299 require_field = RequireField
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar_decorators.py in wrapper(*args, **kwds)
111 # same_field=True might trigger an exception otherwise.
112
--> 113 numfield, new_elems, morphism = number_field_elements_from_algebraics(orig_elems, same_field=True, minimal=True)
114
115 elem_dict = dict(zip(orig_elems, new_elems))
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in number_field_elements_from_algebraics(numbers, minimal, same_field, embedded, prec)
2759 real_numbers = []
2760 for v in numbers:
-> 2761 if v._exact_field().is_complex() and real_case:
2762 # the number comes from a complex algebraic number field
2763 embedded_rt = v.interval_fast(RealIntervalField(prec))
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in _exact_field(self)
4428 if isinstance(sd, (ANRational, ANExtensionElement)):
4429 return sd.generator()
-> 4430 self.exactify()
4431 return self._exact_field()
4432
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8180
8181 if op == 'real':
-> 8182 arg.exactify()
8183 rv = (arg + arg.conjugate()) / 2
8184 rv.exactify()
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8429 left = self._left
8430 right = self._right
-> 8431 left.exactify()
8432 right.exactify()
8433 gen = left._exact_field().union(right._exact_field())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8431 left.exactify()
8432 right.exactify()
-> 8433 gen = left._exact_field().union(right._exact_field())
8434 left_value = gen(left._exact_value())
8435 right_value = gen(right._exact_value())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in union(self, other)
3230 op = QQx(op)
3231 # pari_nf = self._field.pari_nf()
-> 3232 pari_nf = self.pari_field()
3233 factors = list(pari_nf.nffactor(op).lift())[0]
3234 x, y = QQxy.gens()
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in pari_field(self)
3134 if self._pari_field is None:
3135 pari_pol = self._field.pari_polynomial("y")
-> 3136 self._pari_field = pari_pol.nfinit(1)
3137 return self._pari_field
3138
cypari2/auto_gen.pxi in cypari2.gen.Gen_base.nfinit()
KeyboardInterrupt:
Suggestions ?
Emmanuel Charpentier
Jan 30, 2022, 12:22:03 AM1/30/22
to sage-support
FWIW, upgrading to 9.5.rc4 and disabling the use of system’s pari and singular gives the same symptoms. Possible suspect :
charpent@zen-book-flip:~$ sage -pip list | grep pari
cypari2 2.1.2
I’m way out of my depth here…
HTH,
Nils Bruin
Jan 30, 2022, 6:49:22 PM1/30/22
to sage-support
On Saturday, 29 January 2022 at 13:51:14 UTC-8 Emmanuel Charpentier wrote:
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in pari_field(self)3134 if self._pari_field is None: 3135 pari_pol = self._field.pari_polynomial("y") -> 3136 self._pari_field = pari_pol.nfinit(1) 3137 return self._pari_field 3138 cypari2/auto_gen.pxi in cypari2.gen.Gen_base.nfinit() KeyboardInterrupt:
nf_init is a perfectly respectable place to hang. After interruption, can you "%debug" and see what the value of pari_pol is? I'd expect that nfinit determines the ring of integers, which means factoring the discriminant. [QQbar shouldn't need the ring of integers of element arithmetic, but this has been stumbled on before: QQbar often tries to find an "optimized" form of the number field, which for small examples is often quite doable and, if you end up doing a LOT of arithmetic in the same field, is often worth the investment. However, because of factoring of the discriminant, it's fully subexponential in complexity, whereas all the things QQbar needs to do are polynomial time. So we already know that QQbar will do things with the wrong theoretical asymptotic complexity. That can come back and bite you, when reality and asymptotics start behaving similarly (which happens eventually ... asymptotically speakinh)
Emmanuel Charpentier
Jan 31, 2022, 11:59:34 AM1/31/22
to sage-support
Thank you, Nils !
Trying to pinpoint things by running a manual solution step by step, it seems that the stumbling oint is the division on two algebraics. I'll try to make a clear minimal example and report it.
Emmanuel Charpentier
Jan 31, 2022, 6:19:49 PM1/31/22
to sage-support
As advertised, an atempt at a minimal (non-)working example :
# Reproducible minimal example
with seed(0): M = matrix(AA, 3, 3, lambda u,v: AA.random_element())
# Working ring
WR = M.base_ring().algebraic_closure()
# A variable to carry the eigenvalues
l = SR.var("l")
# Vector of unknowns for the eigenvectors
V =vector(list(var("v", n=2))+[SR(1)])
# M.eigenvalues does not return. Get them by hand
# Equation determining the eigenspace
ESE = M-l*identity_matrix(WR, 3)
# Eigenvalues
EVa = ESE.det().roots(ring=WR)
# M.eigenvectors_left() doesn't return.
# Set up a system of equations whose roots define the eigenvectors.
Sys = ESE*V
# Solve for each eigenvalue
# solve([u.subs(l==EVa[0][0]) for u in Sys],list(V[:2])) fails.
# Try to solve manually by substitution
# First equation, first variable.
# solve(Sys[0].subs(l==EVa[0][0]), v0) fails.
# Try to solve by explicit algebraic manipulation
E00 = Sys[0].subs(l==EVa[0][0])
# Isolate v0
V00=E0/E00.coefficient(v0) # fast
At this point, tryoing to proint V00‘s value “never returns”. However, a forst interruption (C-C) prints :
sage: V00
C-c C-c---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in minpoly(self)
4492 try:
-> 4493 return self._minimal_polynomial
4494 except AttributeError:
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__getattr__ (build/cythonized/sage/structure/element.c:4754)()
493 """
--> 494 return self.getattr_from_category(name)
495
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.getattr_from_category (build/cythonized/sage/structure/element.c:4866)()
506 cls = P._abstract_element_class
--> 507 return getattr_from_other_class(self, cls, name)
508
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/cpython/getattr.pyx in sage.cpython.getattr.getattr_from_other_class (build/cythonized/sage/cpython/getattr.c:2633)()
360 dummy_error_message.name = name
--> 361 raise AttributeError(dummy_error_message)
362 attribute = <object>attr
AttributeError: 'AlgebraicNumber' object has no attribute '_minimal_polynomial'
During handling of the above exception, another exception occurred:
KeyboardInterrupt Traceback (most recent call last)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__richcmp__ (build/cythonized/sage/structure/element.c:10418)()
1110 return (<Element>self)._richcmp_(other, op)
1111 else:
-> 1112 return coercion_model.richcmp(self, other, op)
1113
1114 cpdef _richcmp_(left, right, int op):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/coerce.pyx in sage.structure.coerce.CoercionModel.richcmp (build/cythonized/sage/structure/coerce.c:20625)()
1979 assert not (isinstance(x, Element) and
1980 (<Element>x)._parent.get_flag(Parent_richcmp_element_without_coercion))
-> 1981 return PyObject_RichCompare(x, y, op)
1982
1983 # Comparing with coercion didn't work, try something else.
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__richcmp__ (build/cythonized/sage/structure/element.c:10394)()
1108 # an instance of Element. The explicit casts below make
1109 # Cython generate optimized code for this call.
-> 1110 return (<Element>self)._richcmp_(other, op)
1111 else:
1112 return coercion_model.richcmp(self, other, op)
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element._richcmp_ (build/cythonized/sage/structure/element.c:10500)()
1112 return coercion_model.richcmp(self, other, op)
1113
-> 1114 cpdef _richcmp_(left, right, int op):
1115 r"""
1116 Basic default implementation of rich comparisons for elements with
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in _richcmp_(self, other, op)
4897 ci1 = self._value.imag().abs()
4898 ci2 = other._value.imag().abs()
-> 4899 if ci1.overlaps(ci2) and self.minpoly() == other.minpoly():
4900 c = cmp_elements_with_same_minpoly(self, other, self.minpoly())
4901 if c is not None:
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in minpoly(self)
4493 return self._minimal_polynomial
4494 except AttributeError:
-> 4495 self.exactify()
4496 self._minimal_polynomial = self._descr.minpoly()
4497 return self._minimal_polynomial
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8429 left = self._left
8430 right = self._right
-> 8431 left.exactify()
8432 right.exactify()
8433 gen = left._exact_field().union(right._exact_field())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8429 left = self._left
8430 right = self._right
-> 8431 left.exactify()
8432 right.exactify()
8433 gen = left._exact_field().union(right._exact_field())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
4363 if isinstance(od, (ANRational, ANExtensionElement)):
4364 return
-> 4365 self._set_descr(self._descr.exactify())
4366
4367 def _set_descr(self, new_descr):
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in exactify(self)
8431 left.exactify()
8432 right.exactify()
-> 8433 gen = left._exact_field().union(right._exact_field())
8434 left_value = gen(left._exact_value())
8435 right_value = gen(right._exact_value())
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in union(self, other)
3230 op = QQx(op)
3231 # pari_nf = self._field.pari_nf()
-> 3232 pari_nf = self.pari_field()
3233 factors = list(pari_nf.nffactor(op).lift())[0]
3234 x, y = QQxy.gens()
/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py in pari_field(self)
3134 if self._pari_field is None:
3135 pari_pol = self._field.pari_polynomial("y")
-> 3136 self._pari_field = pari_pol.nfinit(1)
3137 return self._pari_field
3138
cypari2/auto_gen.pxi in cypari2.gen.Gen_base.nfinit()
and a second interrupt prints :
KeyboardInterrupt:
Exception ignored in: 'sage.symbolic.expression.py_is_integer'
Traceback (most recent call last):
File "sage/structure/parent.pyx", line 1163, in sage.structure.parent.Parent.__contains__ (build/cythonized/sage/structure/parent.c:10408)
File "sage/structure/parent.pyx", line 898, in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9388)
File "sage/structure/coerce_maps.pyx", line 287, in sage.structure.coerce_maps.NamedConvertMap._call_ (build/cythonized/sage/structure/coerce_maps.c:6085)
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 5011, in _integer_
return AA(self)._integer_(ZZ)
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 5485, in _integer_
self.exactify()
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 4365, in exactify
self._set_descr(self._descr.exactify())
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 8182, in exactify
arg.exactify()
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 4365, in exactify
self._set_descr(self._descr.exactify())
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 8431, in exactify
left.exactify()
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 4365, in exactify
self._set_descr(self._descr.exactify())
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 8431, in exactify
left.exactify()
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 4365, in exactify
self._set_descr(self._descr.exactify())
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 8433, in exactify
gen = left._exact_field().union(right._exact_field())
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 3232, in union
pari_nf = self.pari_field()
File "/usr/local/sage-9/local/var/lib/sage/venv-python3.9/lib/python3.9/site-packages/sage/rings/qqbar.py", line 3136, in pari_field
self._pari_field = pari_pol.nfinit(1)
File "cypari2/auto_gen.pxi", line 23641, in cypari2.gen.Gen_base.nfinit
KeyboardInterrupt:
1.000000000000000?*v0 - 35.57125011095806?*v1 - 492.8896998473554?
Similarly, V00.coefficient(v0).is_one() “never returns”.
The current aritmetic on algebraics is therefore problematic for this kind of problems.
Suggestions ?
Alan Stafford
Feb 2, 2022, 5:57:41 AM2/2/22
to sage-support
What is v0? When I run the above it isn't defined.
Emmanuel Charpentier
Feb 2, 2022, 1:08:20 PM2/2/22
to sage-support
See lines #7-8 of the "Minimal example" :
Le mercredi 2 février 2022 à 11:57:41 UTC+1, alan_thoma...@yahoo.co.uk a écrit :
What is v0? When I run the above it isn't defined.On Monday, January 31, 2022 at 11:19:49 PM UTC Emmanuel Charpentier wrote:As advertised, an atempt at a minimal (non-)working example :
# Reproducible minimal example with seed(0): M = matrix(AA, 3, 3, lambda u,v: AA.random_element()) # Working ring WR = M.base_ring().algebraic_closure() # A variable to carry the eigenvalues l = SR.var("l")
Here is the definition you seek :
# Vector of unknowns for the eigenvectors V =vector(list(var("v", n=2))+[SR(1)])
Is that clearer ?
Alan Stafford
Feb 2, 2022, 2:05:49 PM2/2/22
to sage-support
Do you mean, V =vector(list(var("v0", n=2))+[SR(1)])
Also, what is E0?
Nils Bruin
Feb 2, 2022, 4:15:00 PM2/2/22
to sage-support
On Monday, 31 January 2022 at 15:19:49 UTC-8 Emmanuel Charpentier wrote:
As advertised, an atempt at a minimal (non-)working example :
# Reproducible minimal example with seed(0): M = matrix(AA, 3, 3, lambda u,v: AA.random_element()) # Working ring WR = M.base_ring().algebraic_closure() # A variable to carry the eigenvalues l = SR.var("l") # Vector of unknowns for the eigenvectors V =vector(list(var("v", n=2))+[SR(1)]) # M.eigenvalues does not return. Get them by hand
Actually, for me on 9.5beta9, `M.eigenvalues()` works just fine. So the problem is perhaps just platform-dependent, or there is a very recent change that affected this (my M gets just integer entries from {-2..2})
Looking at the example a bit: you'd be forcing sage to work with a huge compositum if you're actually getting a 3x3 matrix with non-rational algebraic entries: even if they are just independent quadratics, you'd end up in an extension of degree 2^9. This will only work in very limited cases.
One way to get this kind of thing to work is to work with high-precision floats, use numerically (fairly) stable methods to compute the desired answer, and then try to recognize it as algebraic. You probably only care if it is one of fairly low height. You can then try to turn your computation into proof, possibly by tracing through height bounds and showing your precision was sufficient to identify the right solution uniquely.
You could work on automating this kind of thing, but I doubt you'd ever get it to work on a reasonable range of examples; just because the height bounds would be rather ill-behaved.
You can still trace the root cause further on this and perhaps improve arithmetic in AA a bit, but the general shape of the problem you're trying to deal with does not look promising for generally performant methods.
Emmanuel Charpentier
Feb 3, 2022, 1:44:47 AM2/3/22
to sage-support
Le mercredi 2 février 2022 à 22:15:00 UTC+1, Nils Bruin a écrit :
On Monday, 31 January 2022 at 15:19:49 UTC-8 Emmanuel Charpentier wrote:As advertised, an atempt at a minimal (non-)working example :
# Reproducible minimal example with seed(0): M = matrix(AA, 3, 3, lambda u,v: AA.random_element()) # Working ring WR = M.base_ring().algebraic_closure() # A variable to carry the eigenvalues l = SR.var("l") # Vector of unknowns for the eigenvectors V =vector(list(var("v", n=2))+[SR(1)]) # M.eigenvalues does not return. Get them by handActually, for me on 9.5beta9, `M.eigenvalues()` works just fine.
Hmmm… You may have obtained a “less pathological” M than I did, due to possible differences in random numbers generation (notwithstanding my attempt at reproducibility…).
What do you get for M ? I have :
sage: with seed(0): M = matrix(AA, 3, 3, lambda u,v: AA.random_element())
sage: M.apply_map(lambda u:u.radical_expression())
[ -sqrt(2) - 1 -1/4 -2*sqrt(3)]
[ 1/2 1/8*sqrt(33) + 1/8 -1/5*sqrt(29) + 3/5]
[ 0 1/4 1/2]
So the problem is perhaps just platform-dependent, or there is a very recent change that affected this (my M gets just integer entries from {-2..2})
Okay. We have a problem in reproducibility : with seed(0): should entail a reproducible, platform-independent result. It did not. BTW, what is your platform ?
Suggestions on how to document this and file a ticket ?
I agree with the rest of your conclusions, but going to numerical approximations then trying to somehow “recognize” the algebraics they are approximations of somehow denies the whole point of working in QQbar…
Alan Stafford
Feb 4, 2022, 7:20:26 AM2/4/22
to sage-support
On Apple Mac the example above runs on the 9.4 kernel using either the 9.4 or 9.5 interface but not on the 9.5 kernel from either Interface.
Emmanuel Charpentier
Feb 5, 2022, 6:48:47 AM2/5/22
to sage-support
What exactly fails in the example ?
Alan Stafford
Feb 5, 2022, 12:05:42 PM2/5/22
to sage-support
M.eigenvalues() never returns.
Emmanuel Charpentier
Feb 6, 2022, 3:53:51 AM2/6/22
to sage-support
That was my initial complaint... ;-)
Alan Stafford
Feb 6, 2022, 5:52:03 AM2/6/22
to sage-support
My point is that it is the same on Apple Mac, that it isn't related to the recent interface and multiprocessing problems but has just appeared with the 9.5 kernel beta release as it works with the 9.4 kernel.
Dima Pasechnik
Feb 6, 2022, 8:28:08 AM2/6/22
to sage-support
Is it the same M for 9.4 as for 9.5?
--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/0103d0cd-ae31-4e6f-861f-ccbbbda26830n%40googlegroups.com.
Alan Stafford
Feb 6, 2022, 11:35:20 AM2/6/22
to sage-support
No, it works with the 9.4 kernel but not with the 9.5 beta.
Emmanuel Charpentier
Feb 6, 2022, 2:13:02 PM2/6/22
to sage-support
Nins and me got different random matrices : his was composed of integers in (-2..2), and Sage could compute its eigenspaces ; mine was :
sage: with seed(0): M = matrix(AA, 3, 3, lambda u,v: AA.random_element())
sage: M.apply_map(lambda u:u.radical_expression())
[ -sqrt(2) - 1 -1/4 -2*sqrt(3)]
[ 1/2 1/8*sqrt(33) + 1/8 -1/5*sqrt(29) + 3/5]
[ 0 1/4 1/2]
which, BTW, shows a problem about the reproducibility of random elements between Sage versions and/or platforms…
Could you post your random matrix ? And could you try :
M.95 = matrix(AA,3,3,[-sqrt(2) - 1, -1/4, -2*sqrt(3), 1/2, 1/8*sqrt(33) + 1/8, -1/5*sqrt(29) + 3/5, 0, 1/4, 1/2])
ES = M95.eigenspaces_left()
and report if you get a result ?
slelievre
Feb 7, 2022, 12:12:26 AM2/7/22
to sage-support
Reply all
Reply to author
Forward
0 new messages