ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 28 Apr 2021 17:37:06 +0200Issue with inversion of complex symbolic array initialized with numpyhttps://ask.sagemath.org/question/56850/issue-with-inversion-of-complex-symbolic-array-initialized-with-numpy/Hello!
I am trying to invert a symbolic matrix that is initialized with a combination of numpy arrays and I get an error. Below i present a simple code that gives the error.
import numpy as np
a = np.zeros((5,5) , dtype = 'complex')
np.fill_diagonal(a,1)
b = var('x')*a
c = matrix(b)
c.inverse()
Gives the error:
> ECL says: THROW: The catch
> MACSYMA-QUIT is undefined.
What I have noticed up until now is that the issue stems from the data type. When I try casting the numpy array to float before turning it into a matrix, it works. However, my actual code makes use of complex coefficients. I believe it might have to do with how the imaginary part is represented in SageMath in comparison to numpy. ( `j` vs `I` )
When I manually create the symbolic array with `I` for the imaginary part, the `.inverse()` has no issue.Wed, 28 Apr 2021 13:31:54 +0200https://ask.sagemath.org/question/56850/issue-with-inversion-of-complex-symbolic-array-initialized-with-numpy/Answer by slelievre for <p>Hello! </p>
<p>I am trying to invert a symbolic matrix that is initialized with a combination of numpy arrays and I get an error. Below i present a simple code that gives the error.</p>
<pre><code>import numpy as np
a = np.zeros((5,5) , dtype = 'complex')
np.fill_diagonal(a,1)
b = var('x')*a
c = matrix(b)
c.inverse()
</code></pre>
<p>Gives the error:</p>
<blockquote>
<p>ECL says: THROW: The catch
MACSYMA-QUIT is undefined.</p>
</blockquote>
<p>What I have noticed up until now is that the issue stems from the data type. When I try casting the numpy array to float before turning it into a matrix, it works. However, my actual code makes use of complex coefficients. I believe it might have to do with how the imaginary part is represented in SageMath in comparison to numpy. ( <code>j</code> vs <code>I</code> )</p>
<p>When I manually create the symbolic array with <code>I</code> for the imaginary part, the <code>.inverse()</code> has no issue.</p>
https://ask.sagemath.org/question/56850/issue-with-inversion-of-complex-symbolic-array-initialized-with-numpy/?answer=56852#post-id-56852Thanks for reporting.
One alternative approach would be, after creating the array `a`,
instead of multiplying by `x` and then converting to a matrix,
to first convert to a matrix and then multiply by `x`.
Create the array `a`:
sage: import numpy as np
sage: a = np.zeros((5, 5), dtype='complex')
sage: np.fill_diagonal(a, 1)
sage: x = SR.var('x')
sage: b = x*a
sage: a
array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j]])
From the question: multiply array by `x`, then convert to matrix:
sage: b
array([[((1+0j))*x, 0j, 0j, 0j, 0j],
[0j, ((1+0j))*x, 0j, 0j, 0j],
[0j, 0j, ((1+0j))*x, 0j, 0j],
[0j, 0j, 0j, ((1+0j))*x, 0j],
[0j, 0j, 0j, 0j, ((1+0j))*x]], dtype=object)
sage: c = matrix(b)
sage: c
[((1+0j))*x 0j 0j 0j 0j]
[ 0j ((1+0j))*x 0j 0j 0j]
[ 0j 0j ((1+0j))*x 0j 0j]
[ 0j 0j 0j ((1+0j))*x 0j]
[ 0j 0j 0j 0j ((1+0j))*x]
sage: c.parent()
Full MatrixSpace of 5 by 5 dense matrices over Symbolic Ring
sage: c.inverse()
Traceback (most recent call last)
...
TypeError: ECL says: THROW: The catch MACSYMA-QUIT is undefined.
Alternative approach: convert to matrix, then multiply by `x`:
sage: b = matrix(a)
sage: b
[1.0 0.0 0.0 0.0 0.0]
[0.0 1.0 0.0 0.0 0.0]
[0.0 0.0 1.0 0.0 0.0]
[0.0 0.0 0.0 1.0 0.0]
[0.0 0.0 0.0 0.0 1.0]
sage: b.parent()
Full MatrixSpace of 5 by 5 dense matrices over Complex Double Field
sage: c = x * b
sage: c
[1.0*x 0.0 0.0 0.0 0.0]
[ 0.0 1.0*x 0.0 0.0 0.0]
[ 0.0 0.0 1.0*x 0.0 0.0]
[ 0.0 0.0 0.0 1.0*x 0.0]
[ 0.0 0.0 0.0 0.0 1.0*x]
sage: c.parent()
Full MatrixSpace of 5 by 5 dense matrices over Symbolic Ring
sage: c.inverse()
[1.0/x 0 0 0 0]
[ 0 1.0/x 0 0 0]
[ 0 0 1.0/x 0 0]
[ 0 0 0 1.0/x 0]
[ 0 0 0 0 1.0/x]
Wed, 28 Apr 2021 14:45:50 +0200https://ask.sagemath.org/question/56850/issue-with-inversion-of-complex-symbolic-array-initialized-with-numpy/?answer=56852#post-id-56852Comment by yorgos_sot for <p>Thanks for reporting.</p>
<p>One alternative approach would be, after creating the array <code>a</code>,
instead of multiplying by <code>x</code> and then converting to a matrix,
to first convert to a matrix and then multiply by <code>x</code>.</p>
<p>Create the array <code>a</code>:</p>
<pre><code>sage: import numpy as np
sage: a = np.zeros((5, 5), dtype='complex')
sage: np.fill_diagonal(a, 1)
sage: x = SR.var('x')
sage: b = x*a
sage: a
array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j]])
</code></pre>
<p>From the question: multiply array by <code>x</code>, then convert to matrix:</p>
<pre><code>sage: b
array([[((1+0j))*x, 0j, 0j, 0j, 0j],
[0j, ((1+0j))*x, 0j, 0j, 0j],
[0j, 0j, ((1+0j))*x, 0j, 0j],
[0j, 0j, 0j, ((1+0j))*x, 0j],
[0j, 0j, 0j, 0j, ((1+0j))*x]], dtype=object)
sage: c = matrix(b)
sage: c
[((1+0j))*x 0j 0j 0j 0j]
[ 0j ((1+0j))*x 0j 0j 0j]
[ 0j 0j ((1+0j))*x 0j 0j]
[ 0j 0j 0j ((1+0j))*x 0j]
[ 0j 0j 0j 0j ((1+0j))*x]
sage: c.parent()
Full MatrixSpace of 5 by 5 dense matrices over Symbolic Ring
sage: c.inverse()
Traceback (most recent call last)
...
TypeError: ECL says: THROW: The catch MACSYMA-QUIT is undefined.
</code></pre>
<p>Alternative approach: convert to matrix, then multiply by <code>x</code>:</p>
<pre><code>sage: b = matrix(a)
sage: b
[1.0 0.0 0.0 0.0 0.0]
[0.0 1.0 0.0 0.0 0.0]
[0.0 0.0 1.0 0.0 0.0]
[0.0 0.0 0.0 1.0 0.0]
[0.0 0.0 0.0 0.0 1.0]
sage: b.parent()
Full MatrixSpace of 5 by 5 dense matrices over Complex Double Field
sage: c = x * b
sage: c
[1.0*x 0.0 0.0 0.0 0.0]
[ 0.0 1.0*x 0.0 0.0 0.0]
[ 0.0 0.0 1.0*x 0.0 0.0]
[ 0.0 0.0 0.0 1.0*x 0.0]
[ 0.0 0.0 0.0 0.0 1.0*x]
sage: c.parent()
Full MatrixSpace of 5 by 5 dense matrices over Symbolic Ring
sage: c.inverse()
[1.0/x 0 0 0 0]
[ 0 1.0/x 0 0 0]
[ 0 0 1.0/x 0 0]
[ 0 0 0 1.0/x 0]
[ 0 0 0 0 1.0/x]
</code></pre>
https://ask.sagemath.org/question/56850/issue-with-inversion-of-complex-symbolic-array-initialized-with-numpy/?comment=56855#post-id-56855Thank you so much. I would like to also add something that I had to modify in my code.
Because I use python, I had to also define exponential variables as:
import sage.all as sg
e_Iphi = sg.exp(sg.I * sg.var('phi', domain='real') )Wed, 28 Apr 2021 17:37:06 +0200https://ask.sagemath.org/question/56850/issue-with-inversion-of-complex-symbolic-array-initialized-with-numpy/?comment=56855#post-id-56855