# 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!

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.

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, ..... })

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.

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.

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