linsolve bug

[Albert Chan]

Tested on HP-Prime emulator 2.1.14181 (2018 10 16), and XCas 1.9.0-31 (mingw win32)

Example taken from Accuracy of Quadratic Regression, solving quadratic regression coefficients.
Matrix A reduced to a diagonal matrix. We can read off the answer, coefficients = b/a

Cas> a := [1.53963596424e50, 9.47961058243e35, 7.07482042846e19]; /* diagonals */
Cas> b := [-8.42134216805e50,2.08648957898e40,-1.56645830033e27]; /* last column */

Cas> linsolve(diag(a), b)      → []      // ???
Cas> linsolve(diag(a*1), b*1)           // why *1 make a difference?

[−5.4696969697, 22010.2878788, −22141315.3333]

Cas> a := a*1
Cas> b := b*1
Cas> linsolve(diag(a), b)      → []      // *1 not work now ...

[parisse]

This comes from the default value of epsilon: 9e35 is considered to be 0 when compared to 1.5e50 and the matrix is viewed like a singular matrix (2nd row==0).
If you change epsilon to a smaller value (e.g. 1e-30), you will get what you expect.

[Albert Chan]

Thanks for the tip! For my example, smaller epsilon help XCas, but not HP Prime Cas

XCas> epsilon := 1e-30
XCas> linsolve(diag(a), b)

[-5.4696969697, 22010.2878788, -22141315.3333]

Cas> epsilon := 1e-30
Cas> linsolve(diag(a), b)

[C_0, C_1, C_2]

This symbolic result appears if epsilon is 1e-16 or smaller.

[parisse]

The public firmware version of giac is probably too old. It works with an up to date giac version (I checked on the emulator).

[Albert Chan]

Some of linsolve accuracy issue may be my upgrade from XCas 1.5.0 to 1.9.0

Float precision is reduced, from 53-bits mantissa IEEE double, to 48 (what HP Prime Cas has)
Worse, XCas 1.9.0 now round the same as HP Prime Cas too ... truncation!

Cas> h := 2^48 - 1.0
Cas> 2.96+h-h      → 1.03125      // same as XCas 1.9.0
Cas> 2.97+h-h      → 3.03125      // same as XCas 1.9.0

With ill-conditioned matrix, lost of precision and losing round-to-nearest matters!

Code:

/* sx[k] = sum(X.^k), syx[k] = sum(Y.*X.^k) */

sx := float([10,20075,40300645,80903876075,162415528660693]);   
syx := float([11060, 22207030, 44588891682]);

a := [[sx[4]-sx[2]*sx[2]/sx[0], sx[3]-sx[1]*sx[2]/sx[0]], 
      [sx[3]-sx[1]*sx[2]/sx[0], sx[2]-sx[1]*sx[1]/sx[0]]];
b := [syx[2]-syx[0]*sx[2]/sx[0], 
      syx[1]-syx[0]*sx[1]/sx[0]];
epsilon := 1e-30;
linsolve(a,b);

XCas 1.5.0 --> [-5.46969697020, 22010.2878808]
XCas 1.9.0 --> [-5.46452227167, 21989.5114662]

[parisse]

Giac is using 48 bits mantissa by default, unless it is compiled with the -DDOUBLEVAL flag. At least on linux. On windows, I don't remember, perhaps it is different for the 32 bits cygwin version, 64 bits cygwin version and 64 bits mingw version ("native").
This is done in order to have 5 bits to store the type of a gen (generic C type to represent an object) and keep sizeof(gen)=8 bytes=64 bits. If -DDOUBLEVAL is used, then sizeof(gen)=12 bytes on 32 bits architectures and =16 bytes on 64 bits architectures. This means that objects are larger, e.g. matrices take twice more space, which slow down computations (especially when the L1 cache limit is reached).
48 bits is sufficient except for bad conditionned computations, and if you need more then you can calculate with longfloats using e.g. evalf(.,30) to have 30 digits. Multiprecision floats are currently not available on the Prime, they could be sometimes in the future, with quickjs, already available on other calculator ports of giac.

[Albert Chan]

(10-04-2018 07:22 PM)parisse Wrote:  Xcas is using longfloats (via MPFR) for representation of floats having more than 14 digits.

I totally forget about XCas longfloats !

Digits := 16 (or via Setup), and rerun code unchanged, it matched IEEE double accuracy.
BTW, last number is in normalized float notation.

XCas-1.9.0 --> [-5.469696970203201, 0.2201028788082040e5]

Update: MPFR bits = ceil(digits / log10(2)) --> 16 digits == 54 bits.

[Albert Chan]

FYI, this is result of time(linsolve(a,b)). Timing 2nd digit is just a rough guess.

Code:

XCas-1.5.0:         TIME
    53 bits double  6.0E-5
    16 digits MPFR  2.8E-4
    32 digits MPFR  3.2E-4
    
XCas-1.9.0:
    48 bits trunc   6.2E-6
    16 digits MPFR  1.2E-4
    32 digits MPFR  1.3E-4

MPFR timings are basically the same, regardless of setting digits to 16 or 32, or even 100.

Order of magnitude difference between 48 bits and IEEE double is odd.
Perhaps comparing wildly different XCas versions is not a good idea.
But, that's the only 2 versions I had to test.

Anyone has XCas 1.90 (win32), compiled with the -DDOUBLEVAL flag ?

[parisse]

Your systems are too small to be able to deduce something meaningfull. You should try with a randmatrix of order a few hundreds to be able to see something.

[Albert Chan]

100×100 matrix, speed ratio dropped to 2X, which may be due to XCas versions alone

randseed(0);
m := float(ranm(100,100));
time(inverse(m))

Note: both XCas versions produce same random matrix, and can be inverted.

Code:

XCas-1.5.0  TIME
digits:=12  0.011  (53 bits double)
digits:=16  16.    (54 bits MPFR)
digits:=32  18.    (107 bits MPFR)

XCas-1.9.0  TIME
digits:=12  0.0055 (48 bits)
digits:=16  6.4    (54 bits MPFR)
digits:=32  7.4    (107 bits MPFR)

[Albert Chan]

(08-24-2023 02:24 PM)parisse Wrote:  This comes from the default value of epsilon: 9e35 is considered to be 0 when compared to 1.5e50

It may be better *not* to flush numbers to zero. It may make things worse.
Perhaps a warning of possibly rounding issues enough?

(08-25-2023 12:06 AM)Albert Chan Wrote:  \(\left[\begin{array}{cccc}162415528660693&80903876075&40300645&44588891682\\80903876075&40300645&20075&22207030\\40300645&20075&10&11060\end{array}\right]\)

Cas> m := float(m);                 // work with float ... no cheating
Cas> m(1) -= 2007 * m(2);
Cas> m(2) -= 2007 * m(3);

Matrix Column Operations Wrote:We could apply same column operations on top. Continued on from above:

Cas> m := transpose(m)
Cas> m(1) -= 2007 * m(2);
Cas> m(2) -= 2007 * m(3);
Cas> m := transpose(m)

\(\left(\begin{array}{cccc}
342887248.0 & 170720.0 & 10120.0 & 19382472.0 \\
170720.0 & 85.0 & 5.0 & 9610.0 \\
10120.0 & 5.0 & 10.0 & 11060.0
\end{array}\right) \)

Cas> rref(m)             /* WRONG !!! */

\(\left(\begin{array}{cccc}
1.0 & 0.0 & 0.0 & 4.69885643897 \\
0.0 & 1.0 & 0.0 & -9385.88159373 \\
0.0 & 0.0 & 1.0 & 1043.69808063
\end{array}\right) \)

We undo column operations to get back (a,b,c) ...

I was writing this post on column operations, but the answer was wrong!
Column operations does not change values of true a = -361/66 = -5.4(69)

Reduced m, with nice small numbers, seems innocently enough.
Without warning, rref(m) flushed some numbers to zero, and give wrong result.

Making epsilon smaller help ... but it is hard to know ahead of time.

Cas> epsilon := 1e-30      /* was 1e-12 */
Cas> rref(m)

\(\left(\begin{array}{cccc}
1.0 & 0.0 & 0.0 & -5.46969697435 \\
0.0 & 1.0 & 0.0 & 11032.6060699 \\
0.0 & 0.0 & 1.0 & 1125.03030307
\end{array}\right) \)

We undo column operations, to get original variables (a, b, c)

Cas> r := col(Ans, 0)     
Cas> r(3) -= 2007 * r(2)
Cas> r(2) -= 2007 * r(1)

[−5.46969697435, 22010.2878975, −22141315.3521]