Talk:Defactoring algorithms: Difference between revisions

re
No edit summary
 
(2 intermediate revisions by 2 users not shown)
Line 59: Line 59:
from sympy import Matrix, normalforms
from sympy import Matrix, normalforms


def __hnf (main):
def __hnf_col (main):
     return np.flip (np.array (normalforms.hermite_normal_form (Matrix (np.flip (main)).T).T, dtype = int))
     return np.flip (np.array (normalforms.hermite_normal_form (Matrix (np.flip (main))), dtype = int))


def __sat (main):
def __sat (main):
     return np.rint (linalg.inv (__hnf (main.T).T) @ main).astype (int)
     return np.rint (linalg.inv (__hnf_col (main)) @ main).astype (int)
</syntaxhighlight>
</syntaxhighlight>


Line 114: Line 114:


:::::: [[User:FloraC|FloraC]] ([[User talk:FloraC|talk]]) 21:43, 10 February 2023 (UTC)
:::::: [[User:FloraC|FloraC]] ([[User talk:FloraC|talk]]) 21:43, 10 February 2023 (UTC)
::::::: I am actively trying to understand what you're saying. I'm sorry I missed the point you were trying to make with the GitHub thread. Given how long and obscure it was, I should have asked you to be more specific, rather than try to guess at your intentions. It looked to me, based on a quick look over their code and Cohen's work it was based on, that their algorithm worked in a way fairly similar to my existing understanding of how HNF is computed, whereby tracking unimodular transformations is cheap. It's quite possible that their algorithm looks superficially similar but is completely different. I am not interested enough in this problem to figure it out, though, at least not now. Sorry. So again, it's fine by me for now if you want to keep both algorithms.
::::::: But I'm still at a loss about how you've gotten the impression that column Hermite defactoring involves a pseudoinverse step. As you can see in the diagram I shared in my previous reply, it does not. The word "pseudoinverse" only appears once in the article, in a tangential footnote. On my end, "column Hermite defactoring" is a method which Dave Keenan and I developed about a year ago, and the only information about it exists on this wiki article. Are you confusing it with something else? These are the steps:
::::::::# find the transpose A = V<sup>T</sup>;
::::::::# find the unimodular matrix U = First[HermiteDecomposition[A]];
::::::::# find the inverse U<sup>-1</sup>;
::::::::# find the transpose (U<sup>-1</sup>)<sup>T</sup>;
::::::::# take the first ''r'' rows and return.
::::::: You can find an example of us computing the entire thing by hand in the article. No pseudoinverse is required. --[[User:Cmloegcmluin|Cmloegcmluin]] ([[User talk:Cmloegcmluin|talk]]) 22:12, 10 February 2023 (UTC)
:::::::: Oh, allow me to correct myself. After some trial and errors it turns out the pseudoinverse doesn't get you the unimodular matrix, but a singular matrix, so it doesn't work. The correct way seems to be using the augmented matrix to track the elementary row operations in HNF: (A|I) → (H|U), as you pointed out and as one of them showed in the development thread. Here's my implementation of the column Hermite method:
:::::::: <syntaxhighlight lang="python">
def __sat (main):
    r = Matrix (main).rank ()
    unimodular = __hnf_col (np.vstack ((main, np.eye (main.shape[1], dtype = int))))[r:]
    return np.rint (linalg.inv (unimodular)[:r]).astype (int)
</syntaxhighlight>
:::::::: As I said, "Hermite decomposition" isn't a single operation. Here it involves concatenation, HNF, and slice, tho it doesn't involve the pseudoinverse. Compare Pernet-Stein:
:::::::: <syntaxhighlight lang="python">
def __sat (main):
    r = Matrix (main).rank ()
    return np.rint (linalg.inv (__hnf_col (main)[:, :r]) @ main).astype (int)
</syntaxhighlight>
:::::::: Note the slice and thus computing the rank are redundant; they're there just for futureproofness.
:::::::: [[User:FloraC|FloraC]] ([[User talk:FloraC|talk]]) 10:58, 11 February 2023 (UTC)
:::::::: Edit: update the code to adopt column-style HNF so as to get rid of all the transpositions.
:::::::: [[User:FloraC|FloraC]] ([[User talk:FloraC|talk]]) 07:08, 15 February 2023 (UTC)
Return to "Defactoring algorithms" page.