Add implementation explanation in computer_gcd

Author: PierreQuinton

src/torchjd/sparse/_linalg.py

@@ -135,6 +135,16 @@ def compute_gcd(S1: Tensor, S2: Tensor) -> tuple[Tensor, Tensor, Tensor]:
     S2 = G @ K2
     with G having minimal rank.
 
+    Implementation logic:
+    The concatenated matrix [S1 | S2] spans exactly the sum of the lattices generated by S1 and S2.
+    This is because S1 @ u1 + S2 @ u2 = [S1 | S2] @ [u1.T | u2.T].T
+    The reduced HNF decomposition of [S1 | S2] yields G, U, V where the G.shape[1] is the rank of
+    [S1 | S2] and [S1 | S2] = G @ V. This means that
+    S1 = G @ V[:, :m1]
+    S2 = G @ V[:, m1:]
+    This is the target common factorization. It is the greatest as the lattice spanned by G is the
+    same as that spanned by [S1 | S2].
+
     Args:
         S1, S2: torch.Tensors (m x n1), (m x n2)
 
@@ -149,22 +159,6 @@ def compute_gcd(S1: Tensor, S2: Tensor) -> tuple[Tensor, Tensor, Tensor]:
     A = torch.cat([S1, S2], dim=1)
     G, U, V = hnf_decomposition(A)
 
-    # H = [S1 | S2] @ U
-    # [S1 | S2] = H @ V
-    #
-    # S1 = H @ V[:, :m1]
-    # S2 = H @ V[:, m1:]
-    #
-    # K1 = V[:, :m1]
-    # K2 = V[:, m1:]
-    # G = H
-    #
-    # S1 = G @ K1
-    # S2 = G @ K2
-    #
-    # SLT(p1, S1) = SLT(SLT(p1, K1), G)
-    # SLT(p2, S2) = SLT(SLT(p2, K2), G)
-
     K1 = V[:, :n1]
     K2 = V[:, n1:]