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 <[email protected]> wrote: > On Tue, 11 Jan 2022 at 12:30, Gerardo Suarez <[email protected]> > 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 [email protected]. > 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CABKqA0ZOY0hRG%2B1XsTeSiE0UZK9O88DYnc2rzeM4vWNBQwxB-A%40mail.gmail.com.
