matrix is found singular over CC, nonsingular over CDF

I am trying to take the inverse of complex valued matrix where every entry is either +/-1,+/-i,or 0. The matrices are quite large and I am getting substantial rounding errors when I construct the matrix using "matrix(CDF,array)" and then take the inverse, but I still get an output.

When I construct the matrix using the same array but using "matrix(CC,array)" and then take the inverse I get this error: ZeroDivisionError: input matrix must be nonsingular

I am confused because the array when I use "matrix(CC,array)" is exactly the same as the array when I use CDF, however "matrix(CDF,array)" is nonsingular and the function inverse works. Why is this happening and is there another way I can avoid large rounding errors when taking the inverse?

I hope this is enough detail. I'm not really sure how to provide an example of this error occurring since the smallest example of the array I am working with is still quite large (52x52) and I haven't been able to get this error to occur with smaller simpler arrays.

This is the encoded matrix over the Gaussian Integers - I had been using QQbar() to run it since it worked to get the inverse.

b = [286102294921876, 178813934326175001, 69849193096160908206250, 27284841053187847167978515625, 10658141036401502788066894531250000, 179100036621093756, 111937522888183596880, 69960951805114746112500, 43725594878196716328140625, 27343048714101314544736328125, 17080310499295592308075000000000, 10675194062059745192527954101562500, 1716613769531255, 112652778625488296880, 44004991650581359960953125, 17189449863508343696746875000000, 63948846218409016728401336669921875, 895500183105468780, 559687614440917984400, 349804759025573730562500, 218627974390983581640703125, 136715243570506572723681640625, 85401552496477961540375000000000, 53375970310298725962639770507812500, 8583068847656275, 563263893127441484400, 85947249317541718483734375000000, 319744231092045083642006683349609375, 4477500915527343900, 2798438072204589922000, 1749023795127868652812500, 1093139871954917907959375000, 683212419971823692323242187500, 427007762482389807701875000000000, 266879851551493629813198852539062500, 42915344238281375, 2816319465637207422000, 1100124791264533997803125000, 429736246587708592418671875000000, 1598721155460225418210033416748046875, 22387504577636719500, 13992190361022949610000, 8745118975639343264062500, 5465699359774589539796875000, 3416062099859118461616210937500, 2135038812411949038509375000000000, 1334399257757468149065994262695312500, 35762786865235000, 13969838619232179688125, 5456968210637569433595703125, 2131628207280300557632452392578125, 1332267629550187848508377075195312500]

The issue I'm now having is that I expected the inverse this matrix to look much different - as well as the determinant to be much larger.

edit retag close merge delete

Thanks for adding the encoded example. It seems to reveal several bugs.

What remains to do now:

• find a minimal example to illustrate each of the bugs
• check whether these bugs are already tracked on the Sage Trac server
• if not, file them there
• and of course solve them!
( 2021-05-18 19:03:59 +0200 )edit

What are the bugs? I cannot reproduce the ZeroDivisionError in Sage 9.3. The matrix is well-conditioned, so I would expect inversion to work without much problem.

( 2021-05-20 21:27:11 +0200 )edit

Note that your matrix is quite sparse, so you can use .dict() to write it more directly as

matrix(GaussianIntegers(), {(0, 0): 1, (0, 20): I, (1, 0): 1, (1, 5): 1, ..... })
( 2021-05-20 21:33:39 +0200 )edit

Good point about the dict approach since the matrix is quite sparse.

I also cannot reproduce the ZeroDivisionError but by exploring the example in my answer I isolate three things I think are bugs.

Certainly the fact that, after changing ring to CC, the matrix has determinant 512.0*I instead of 2 is a bug.

( 2021-05-21 01:30:23 +0200 )edit

Oh, I had missed that part of the answer. Indeed, the incorrect determinant over CC looks quite bad (in most cases, CDF is a better choice), but I think this boils down to numerical issues. The determinant is computed via ._charpoly_df(), the division-free implementation. I am not familiar with the algorithm, but I expect it computes products with 52 factors. At least, the computation attains intermediate values larger than 1e19, which is a problem if the precision is not increased.

Alternatively, .determinant(algorithm='hessenberg') gives correct results (if the determinant has not already been cached). Actually, I have observed before that the hessenberg implementation can be much faster as well.

( 2021-05-21 20:23:38 +0200 )edit

Sort by ยป oldest newest most voted

Dealing with rounding errors

Why is this happening and is there another way I can avoid large rounding errors when taking the inverse?

This is happening because CC and CDF use floating-point numbers, and computations are thus inexact and subject to rounding errors.

Instead of working with CC or CDF, one could

• use matrices over GaussianIntegers() or ZZ[i] or QQ[i] or QQbar for exact computations

• use matrices over CBF ("complex ball field") for floating-point computations with rigorous error bounds

Providing a large example

I'm not really sure how to provide an example of this error occurring since the smallest example of the array I am working with is (52x52).

Communicating a 52 x 52 matrix with entries in {0, 1, -1, i, -i} is doable.

You can encode the matrix and provide decoding instructions rather than the matrix itself.

One possible encoding, since the entries can only take 5 values, take them as 5-ary digits; then each row can be thought of as the 5-ary expansion of a number. If we agree on the digits and communicate the numbers, we're good.

Here is a random matrix a with entries in {0, 1, -1, i, -i}:

sage: G = GaussianIntegers()
sage: U = [0, 1, -1, i, -i]

sage: a = matrix(G, 52, lambda i, j: U[randint(0, 4)])


Encoded matrix:

sage: b = [Integer(tuple(U.index(e) for e in row), 5) for row in a]

sage: print(b)
[1986835196674783292709127771671427698, 372532889760445898415871372893604877,
29903254311383361818758881193795318, 1358695302115484662298806740046547384,
1537321601024823802437657316988965568, 807200500103255107975400681439180354,
567350510853383980331713019521919561, 1502850803304703506587108729499355006,
846884730762028905888265726139310683, 137381915837271267181515478275643286,
522105367460276163750943030646864364, 1062521690576328608603843582516522599,
899547582492594668160772320652245190, 767732247768012897653127687752209191,
629354811960707798077793065391805161, 922202162939977177914455168534168765,
1729869497057901779087386639626774141, 1868018985439944039499425518652161649,
1296299082686624643142119402841062894, 1064799841987832209268578137251968874,
2155066789175869415783929309121456698, 267366791364966820259349410583047599,
1632002494019334195780113340724191796, 2192661273536615612340019587291366905,
1762917404578120599684457836032844498, 228986171419406644876547615831863389,
2101920423365829833128536344495944345, 660985714605962040665526677335410210,
873334831095240711598788236281716427, 1744121894438415469560277438117935725,
2210059965689974187926035573546483126, 1594508474721458635473949546648185838,
518211385383214049244941993118895991, 231206763245123017363282654207827645,
1748059359372987512622900717899928645, 259988738192311572545249238464600548,
2007721547194093983826996655201647232, 1240843285105582949996674635390136465,
1736982326629685244276013105887136315, 435682393610221223713098672542297123,
66566788503415822985489945441756435, 1141561970732274264554486688206105649,
2044822779971527183978088623745613537, 17299608296692077102235770236630252,
780853597789891512314851615527578148, 427428566224687826328803693567479441,
1922351158930802787176924223394683886, 2104373328835002410314401979576425296,
1317300586469460326137557530572476230, 1736284365530465805899187629996562866,
2132757211837100962505268985637520292, 147031314086172436786678135727403627]


To decode:

sage: aa = matrix(G, 52, [[U[d] for d in r.digits(5, padto=52)] for r in b])
sage: aa == a
True


Here is a function that takes a matrix as an argument, and returns "short enough code" to recreate that matrix.

def matrix_code(a):
r"""
Return code for a matrix with entries in {0, 1, -1, i, -i}.

EXAMPLE::

sage: G = GaussianIntegers()
sage: U = [0, 1, -1, i, -i]
sage: a = matrix(G, 52, lambda i, j: U[randint(0, 4)])
sage: matrix_code(a)
...
"""
G = GaussianIntegers()
U = [0, 1, -1, i, -i]
m = a.nrows()
n = a.ncols()
b = [Integer(tuple(U.index(e) for e in row), 5) for row in a]
print(f"""
G = GaussianIntegers()
U = [0, 1, -1, i, -i]
b = {b}
a = matrix(G, {m}, {n}, [[U[d] for d in r.digits(5, padto={n})] for r in b])
""")


Run it on your matrix. Edit your question, pasting the output. Then we can all play with the same matrix.

Exploring an example

(Added after the question was updated to provide a concrete example.)

Thanks for the concrete example. It seems to uncover several bugs!

Define the gaussian integers and their fraction field:

sage: Zi = GaussianIntegers()
sage: Qi = Zi.fraction_field()


The "digits" for decoding the encoded matrix in the updated question:

sage: U = [0, 1, -1, i, -i]


The decoded matrix:

sage: a = matrix(Zi, 52, [[U[d] for d in r.digits(5, padto=52)] for r in b])


Several variations of the same matrix, varying the base ring:

sage: a_Zi = a
sage: a_Qi = a.change_ring(Qi)
sage: a_CC = a.change_ring(CC)
sage: a_CDF = a.change_ring(CDF)
sage: a_QQbar = a_Qi.change_ring(QQbar)


First bug: changing ring directly from the gaussian integers to QQbar fails:

sage: a.change_ring(QQbar)
Traceback (most recent call last)
...
AttributeError: 'AbsoluteOrder_with_category' object has no attribute 'polynomial'


Second bug: all our matrices don't compare equal:

sage: aa = [a_Zi, a_Qi, a_QQbar, a_CC, a_CDF]
sage: print('\n'.join(' '.join(str(int(a == b)) for b in aa) for a in aa))
1 1 0 1 1
1 1 0 1 1
0 0 1 1 1
1 1 1 1 1
1 1 1 1 1


Third bug: the determinant computation over CC fails:

sage: ringdets = [('ZZ[i]', a_Zi.det()),
....:             ('QQ[i]', a_Qi.det()),
....:             ('QQbar', a_QQbar.det()),
....:             ('CC', a_CC.det()),
....:             ('CDF', a_CDF.det())]

sage: print('\n'.join('%5s %18s' % ringdet for ringdet in ringdets))
ZZ[i]                  2
QQ[i]                  2
QQbar                  2
CC 512.000000000000*I
CDF                2.0


Taking the inverse

I don't observe problems when taking the inverse:

sage: b_Zi = a_Zi.inverse()
sage: b_Qi = a_Qi.inverse()
sage: b_QQbar = a_QQbar.inverse()
sage: b_CDF = a_CDF.inverse()
sage: b_CC = a_CC.inverse()


and the inverse seems to be computed correctly; the result for CC is rather more accurate than for CDF.

sage: print(b_Zi[:8, :4])
[  -1    1    0    0]
[  -2    1   -1   -1]
[   0    0   -1   -1]
[   0    0    0  1/2]
[   0    0    0  1/2]
[   0    1   -1   -1]
[  -1    1   -1   -1]
[   1    0   -1 -1/2]
sage: print(b_Qi[:8, :4])
[  -1    1    0    0]
[  -2    1   -1   -1]
[   0    0   -1   -1]
[   0    0    0  1/2]
[   0    0    0  1/2]
[   0    1   -1   -1]
[  -1    1   -1   -1]
[   1    0   -1 -1/2]
sage: print(b_QQbar[:8, :4])
[  -1    1    0    0]
[  -2    1   -1   -1]
[   0    0   -1   -1]
[   0    0    0  1/2]
[   0    0    0  1/2]
[   0    1   -1   -1]
[  -1    1   -1   -1]
[   1    0   -1 -1/2]
sage: print(b_CDF[:8, :4])
[    -0.9999999999999998      1.0000000000000004  -7.771561172376096e-16  -2.220446049250313e-16]
[    -1.9999999999999996      1.0000000000000009                    -1.0     -1.0000000000000004]
[ 1.1102230246251565e-16   5.551115123125783e-17                    -1.0     -0.9999999999999999]
[ -6.106226635438361e-16  2.7755575615628914e-16  2.7755575615628914e-16      0.5000000000000002]
[-1.1102230246251568e-16 -1.1102230246251563e-16                     0.0      0.4999999999999999]
[  5.551115123125783e-16      0.9999999999999997     -1.0000000000000004     -1.0000000000000004]
[    -1.0000000000000004      1.0000000000000009     -1.0000000000000004                    -1.0]
[     0.9999999999999998  1.1102230246251565e-16     -1.0000000000000004     -0.5000000000000002]
sage: print(b_CC[:8, :4])
[ -1.00000000000000   1.00000000000000  0.000000000000000  0.000000000000000]
[ -2.00000000000000   1.00000000000000  -1.00000000000000  -1.00000000000000]
[ 0.000000000000000  0.000000000000000  -1.00000000000000  -1.00000000000000]
[ 0.000000000000000  0.000000000000000  0.000000000000000  0.500000000000000]
[ 0.000000000000000  0.000000000000000  0.000000000000000  0.500000000000000]
[ 0.000000000000000   1.00000000000000  -1.00000000000000  -1.00000000000000]

more