Would Mathematic return a Penrose inverse, if det(M) =0? I do not know
Mathematica at all.
On Tue 11. Jan 2022 at 20:08, Oscar Benjamin <oscar.j.benja...@gmail.com>
wrote:
> On Tue, 11 Jan 2022 at 12:30, Gerardo Suarez <gerardofingu...@gmail.com>
> wrote:
> >
> > Sure, sorry about the image, I was unaware of how easy is to generate
> python code with the print methods
> >
> > from sympy import Symbol, Matrix,I
> > gamma_c = Symbol('gamma_c')
> > n_c = Symbol('n_c')
> > gamma_h = Symbol('gamma_h')
> > n_h = Symbol('n_h')
> > epsilon_c = Symbol('epsilon_c')
> > g = Symbol('g')
> > epsilon_h = Symbol('epsilon_h')
> > e = Matrix([[-gamma_c*(n_c + 1) - gamma_h*(n_h + 1), 0, 0, 0, 0,
> gamma_c*n_c, 0, 0, 0, 0, gamma_h*n_h, 0, 0, 0, 0, 0], [0, -I*epsilon_c -
> gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*(n_h + 1), I*g, 0, 0, 0, 0,
> 0, 0, 0, 0, gamma_h*n_h, 0, 0, 0, 0], [0, I*g, -I*epsilon_h - gamma_c*(n_c
> + 1) - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, 0, 0, 0, 0, gamma_c*n_c, 0, 0,
> 0, 0, 0, 0, 0, 0], [0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 -
> gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 - I*(epsilon_c + epsilon_h), 0, 0, 0,
> 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, -gamma_c*n_c/2 - gamma_c*(n_c +
> 1)/2 - gamma_h*(n_h + 1) - I*(epsilon_h - (epsilon_c + epsilon_h)), 0, 0,
> 0, -I*g, 0, 0, 0, 0, 0, gamma_h*n_h, 0], [gamma_c*(n_c + 1) - gamma_h*n_h,
> 0, 0, 0, 0, -gamma_c*n_c - gamma_h*n_h - gamma_h*(n_h + 1), I*g, 0, 0,
> -I*g, -gamma_h*n_h, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, I*g, -gamma_c*n_c/2 -
> gamma_c*(n_c + 1)/2 - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 - I*(-epsilon_c +
> epsilon_h), 0, 0, 0, -I*g, 0, 0, 0, 0, 0], [0, 0, gamma_c*(n_c + 1), 0, 0,
> 0, 0, -I*epsilon_h - gamma_c*n_c - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2, 0,
> 0, 0, -I*g, 0, 0, 0, 0], [0, 0, 0, 0, -I*g, 0, 0, 0, -gamma_c*(n_c + 1) -
> gamma_h*n_h/2 - gamma_h*(n_h + 1)/2 - I*(epsilon_c - (epsilon_c +
> epsilon_h)), 0, 0, 0, 0, gamma_c*n_c, 0, 0], [0, 0, 0, 0, 0, -I*g, 0, 0, 0,
> -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2
> - I*(epsilon_c - epsilon_h), I*g, 0, 0, 0, 0, 0], [-gamma_c*n_c +
> gamma_h*(n_h + 1), 0, 0, 0, 0, -gamma_c*n_c, -I*g, 0, 0, I*g, -gamma_c*n_c
> - gamma_c*(n_c + 1) - gamma_h*n_h, 0, 0, 0, 0, 0], [0, gamma_h*(n_h + 1),
> 0, 0, 0, 0, 0, -I*g, 0, 0, 0, -I*epsilon_c - gamma_c*n_c/2 - gamma_c*(n_c +
> 1)/2 - gamma_h*n_h, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
> -gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 - gamma_h*n_h/2 - gamma_h*(n_h + 1)/2
> + I*(epsilon_c + epsilon_h), 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0,
> gamma_c*(n_c + 1), 0, 0, 0, 0, I*epsilon_h - gamma_c*n_c - gamma_h*n_h/2 -
> gamma_h*(n_h + 1)/2, I*g, 0], [0, 0, 0, 0, gamma_h*(n_h + 1), 0, 0, 0, 0,
> 0, 0, 0, 0, I*g, I*epsilon_c - gamma_c*n_c/2 - gamma_c*(n_c + 1)/2 -
> gamma_h*n_h, 0], [gamma_c*n_c + gamma_h*n_h, 0, 0, 0, 0, gamma_c*n_c +
> gamma_h*n_h + gamma_h*(n_h + 1), 0, 0, 0, 0, gamma_c*n_c + gamma_c*(n_c +
> 1) + gamma_h*n_h, 0, 0, 0, 0, 0]])
> >
> > However, I care more about how to speed up matrix inversion in general,
> as it was terribly fast in Mathematica and it never finished in sympy. I
> guess there is a more efficient way than using
> >
> > e.inv()
>
> The determinant of this matrix is zero so it isn't invertible. I
> checked this using DomainMatrix:
>
> In [87]: from sympy.polys.matrices import DomainMatrix
>
> In [88]: dM = DomainMatrix.from_Matrix(e)
>
> In [89]: %time ok = dM.charpoly()
> CPU times: user 2min 36s, sys: 100 ms, total: 2min 36s
> Wall time: 2min 36s
>
> In [90]: ok[-1] # constant term of charpoly is the determinant
> Out[90]: (0 + 0*I)
>
> You can check more quickly just by substituting random values for the
> symbols:
>
> In [91]: vals = dict(zip(e.free_symbols, randMatrix(7, 1)))
>
> In [92]: e.subs(vals).det()
> Out[92]: 0
>
> In [93]: vals
> Out[93]: {ε_c: 6, εₕ: 92, g: 44, γ_c: 84, γₕ: 63, n_c: 37, nₕ: 42}
>
> Are you sure that this matrix is the same as the one that you used in
> Mathematica?
>
> Did Mathematica actually return an inverse or did it just say that the
> matrix is not invertible?
>
> Also I can see now how inv can be made much faster for cases like
> this. The matrix has a number of different components:
>
> In [101]: e.connected_components()
> Out[101]: [[0, 5, 6, 10, 9, 15], [1, 2, 7, 11], [3], [4, 8, 13, 14], [12]]
>
> The inverse can be computed from the submatrices generated by the
> numbered rows and columns. The determinant is the product of their
> determinants and the first submatrix has a zero determinant:
>
> In [110]: ccs = e.connected_components()
>
> In [111]: e[ccs[0],ccs[0]].det()
> Out[111]: 0
>
> In fact this is zero because its right column is zero which in turn is
> because column 15 of the matrix is zero:
>
> In [121]: list(e[:,15])
> Out[121]: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
>
> Not sure why I didn't notice that straight away!
>
> In [127]: e.print_nonzero()
> [X X X ]
> [ XX X ]
> [ XX X ]
> [ X ]
> [ X X X ]
> [X XX XX ]
> [ XX X ]
> [ X X X ]
> [ X X X ]
> [ X XX ]
> [X XX XX ]
> [ X X X ]
> [ X ]
> [ X XX ]
> [ X XX ]
> [X X X ]
>
> That's definitely a case that could be handled more efficiently by the
> various routines but it makes me wonder if you showed the right
> matrix.
>
> --
> Oscar
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/CAHVvXxTtngh4LHPu2aid9bqycn%2B7cvg0cZ8Ey0_r5ms4CYzQig%40mail.gmail.com
> .
>
--
Best regards,
Peter Stahlecker
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to sympy+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/CABKqA0ZOY0hRG%2B1XsTeSiE0UZK9O88DYnc2rzeM4vWNBQwxB-A%40mail.gmail.com.
