[sage-devel] Re: Determinant of a SR matrix substituted by a floating number is not numerical stable

2023-01-04 Thread 陳孜穎
Sorry that I missed some important message. The order if save and load are 
not inversed. What I want to do is to load the data stored in the first run 
when running the second time of my program. The intention is I guessed (may 
not be true) some implicit initial states of each run affect the result, 
and I also want to detect the problem caused by the initial states.

John H Palmieri 在 2023年1月5日 星期四凌晨2:53:37 [UTC+8] 的信中寫道:

> Your current code has fragments like this:
>
> weighted_load = load_obj(f'./sym.sobj')
> print(f'weighted == weighted_load {weighted == weighted_load}')
> save(weighted, f'./sym.sobj')
>
> In particular, "sym.sobj" is read before it is written, so 'weighted_load' 
> will be None. If I put the "save" lines earlier, I see this:
>
> weighted == weighted_load True
> realm == realm_load: True
> realm != realm_load: False
> determinant
> -7.116356115499656e-14
> -7.116356115499656e-14
> -7.116356115499656e-14
>
> Having said all of that, I think that determinants tend to be numerically 
> unstable, and it would not be surprising for a number so close to zero to 
> display that instability. A small change in one of the entries of the 
> matrix, due to how the decimals are represented, could easily change the 
> determinant.
>
> On Wednesday, January 4, 2023 at 9:12:55 AM UTC-8 ivana...@gmail.com 
> wrote:
>
>> I found I got different determinant value every time I ran my program. I 
>> tried to save the current matrix and load the matrix in the next run, and I 
>> compared two matrices. The compare results are not stable either. Sometime 
>> they are the same and not the same at the same time (A == B is False and A 
>> != B is False either)!
>>
>> Here is the minimal case I could figure out. I tried to replace the 
>> testing matrix to a simpler one, but the issue could not be reproduced 
>> after replacing. I am sorry about giving you some many lines of code. I 
>> guess it's not easy to read code on Google Groups, so the file is attached.
>>
>> from functools import cache
>>
>> from sage.all_cmdline import *
>>
>>
>> l = 1
>>
>>
>> def load_obj(filename):
>> try:
>> return load(filename)
>> except FileNotFoundError:
>> return None
>>
>>
>> @cache
>> def expected_distance(n, l):
>> if n <= -2:
>> return 0
>> if n == -1:
>> return -2 * x
>> if n == 0:
>> return 1
>> return ((
>> - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
>> - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n 
>> - 2, l)
>> ) / (8 * (n + 1) * x)).full_simplify()
>>
>> # Generate the matrix for testing. I tried to figure out a minimal 
>> reproducible
>> # code, but I failed.
>> depth = 8
>> functions = [expected_distance(i, l) for i in range(depth * 2 - 1)]
>> hankel = matrix.hankel(functions[:depth], functions[depth:depth * 2 - 1], 
>> SR)
>> weighted = matrix(depth, depth, [(hankel[i][j] / hankel[0][j] / hankel[i
>> ][0])
>>  for i in range(depth)
>>  for j in range(depth)])
>>
>> # Only with this useless line, you can see the first and the second cases
>> # listed in the comment (*) below.
>> weighted = weighted.submatrix(0, 0, depth, depth)
>>
>> weighted_load = load_obj(f'./sym.sobj')
>> print(f'weighted == weighted_load {weighted == weighted_load}')
>> save(weighted, f'./sym.sobj')
>>
>> realm = weighted.substitute({x: -0.0070218})
>> realm2 = weighted.substitute({x: -0.0070218})
>> realm_load = load_obj(f'./real.sobj')
>> save(realm, f'./real.sobj')
>>
>> # (*) There are three possible output of this two lines
>> # 1.
>> # realm == realm_load: False
>> # realm != realm_load: False
>> # 2.
>> # realm == realm_load: True
>> # realm != realm_load: False
>> # 3.
>> # realm == realm_load: False
>> # realm != realm_load: True
>> #
>> # The possibilities of each case are almost the same.
>> print(f'realm == realm_load: {realm == realm_load}')
>> print(f'realm != realm_load: {realm != realm_load}')
>>
>> # The determinant values are different in defferent runs
>> print('determinant')
>> print(realm.det())
>> print(realm2.det())
>> print(realm_load.det())
>>
>> # Only print something when the third case happened.
>> for i in range(realm.nrows()):
>> for j in range(realm.ncols()):
>> if realm[i][j] != realm_load[i][j]:
>> print(f'{(i, j)}')
>> print(realm[i][j])
>> print(realm_load[i][j])
>>
>>

-- 
You received this message because you are subscribed to the Google Groups 
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sage-devel+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/sage-devel/52848b5f-a881-4aaa-b603-2fe43224bb29n%40googlegroups.com.


[sage-devel] Re: Determinant of a SR matrix substituted by a floating number is not numerical stable

2023-01-04 Thread John H Palmieri
Your current code has fragments like this:

weighted_load = load_obj(f'./sym.sobj')
print(f'weighted == weighted_load {weighted == weighted_load}')
save(weighted, f'./sym.sobj')

In particular, "sym.sobj" is read before it is written, so 'weighted_load' 
will be None. If I put the "save" lines earlier, I see this:

weighted == weighted_load True
realm == realm_load: True
realm != realm_load: False
determinant
-7.116356115499656e-14
-7.116356115499656e-14
-7.116356115499656e-14

Having said all of that, I think that determinants tend to be numerically 
unstable, and it would not be surprising for a number so close to zero to 
display that instability. A small change in one of the entries of the 
matrix, due to how the decimals are represented, could easily change the 
determinant.

On Wednesday, January 4, 2023 at 9:12:55 AM UTC-8 ivana...@gmail.com wrote:

> I found I got different determinant value every time I ran my program. I 
> tried to save the current matrix and load the matrix in the next run, and I 
> compared two matrices. The compare results are not stable either. Sometime 
> they are the same and not the same at the same time (A == B is False and A 
> != B is False either)!
>
> Here is the minimal case I could figure out. I tried to replace the 
> testing matrix to a simpler one, but the issue could not be reproduced 
> after replacing. I am sorry about giving you some many lines of code. I 
> guess it's not easy to read code on Google Groups, so the file is attached.
>
> from functools import cache
>
> from sage.all_cmdline import *
>
>
> l = 1
>
>
> def load_obj(filename):
> try:
> return load(filename)
> except FileNotFoundError:
> return None
>
>
> @cache
> def expected_distance(n, l):
> if n <= -2:
> return 0
> if n == -1:
> return -2 * x
> if n == 0:
> return 1
> return ((
> - 4 * (2 * (n + 1) - 1) * expected_distance(n - 1, l) \
> - n * ((n + 1) * (n - 1) - 4 * l * (l + 1)) * expected_distance(n 
> - 2, l)
> ) / (8 * (n + 1) * x)).full_simplify()
>
> # Generate the matrix for testing. I tried to figure out a minimal 
> reproducible
> # code, but I failed.
> depth = 8
> functions = [expected_distance(i, l) for i in range(depth * 2 - 1)]
> hankel = matrix.hankel(functions[:depth], functions[depth:depth * 2 - 1], 
> SR)
> weighted = matrix(depth, depth, [(hankel[i][j] / hankel[0][j] / hankel[i][
> 0])
>  for i in range(depth)
>  for j in range(depth)])
>
> # Only with this useless line, you can see the first and the second cases
> # listed in the comment (*) below.
> weighted = weighted.submatrix(0, 0, depth, depth)
>
> weighted_load = load_obj(f'./sym.sobj')
> print(f'weighted == weighted_load {weighted == weighted_load}')
> save(weighted, f'./sym.sobj')
>
> realm = weighted.substitute({x: -0.0070218})
> realm2 = weighted.substitute({x: -0.0070218})
> realm_load = load_obj(f'./real.sobj')
> save(realm, f'./real.sobj')
>
> # (*) There are three possible output of this two lines
> # 1.
> # realm == realm_load: False
> # realm != realm_load: False
> # 2.
> # realm == realm_load: True
> # realm != realm_load: False
> # 3.
> # realm == realm_load: False
> # realm != realm_load: True
> #
> # The possibilities of each case are almost the same.
> print(f'realm == realm_load: {realm == realm_load}')
> print(f'realm != realm_load: {realm != realm_load}')
>
> # The determinant values are different in defferent runs
> print('determinant')
> print(realm.det())
> print(realm2.det())
> print(realm_load.det())
>
> # Only print something when the third case happened.
> for i in range(realm.nrows()):
> for j in range(realm.ncols()):
> if realm[i][j] != realm_load[i][j]:
> print(f'{(i, j)}')
> print(realm[i][j])
> print(realm_load[i][j])
>
>

-- 
You received this message because you are subscribed to the Google Groups 
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to sage-devel+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/sage-devel/a7e473a2-7dfa-4e69-9f15-651ef6e77c09n%40googlegroups.com.