Fix/normalized hypergraph laplacian

tests/linalg/test_matrix.py

@@ -1,9 +1,11 @@
 import numpy as np
 import pytest
+from scipy.linalg import eigh, eigvalsh
 from scipy.sparse import csr_array
 from scipy.sparse.linalg import norm as spnorm
 
 import xgi
+from xgi.exception import XGIError
 
 
 def test_incidence_matrix(edgelist1, edgelist3, edgelist4):
@@ -605,25 +607,57 @@ def test_normalized_hypergraph_laplacian():
 
     assert isinstance(L2, np.ndarray)
     assert np.all(L1.toarray() == L2)
-    assert np.all(np.diag(L2) == 0.5)
 
-    L3, d = xgi.normalized_hypergraph_laplacian(H, index=True)
+    evals = eigh(L2, eigvals_only=True)
+    negative_evals = list(filter(lambda e: e < 0, evals))
+    assert (not negative_evals) or (np.allclose(negative_evals, 0))
 
+    L3, d = xgi.normalized_hypergraph_laplacian(H, index=True)
     assert d == {0: 1, 1: 2, 2: 3, 3: 4, 4: 5, 5: 6, 6: 7, 7: 8}
+
+    # sqrt(d) is eigenvector with eigenvalue 0
+    sqrt_d = np.array([np.sqrt(d) for d in H.nodes.degree.aslist()])
+    assert np.allclose(L3 @ sqrt_d, 0)
+
     true_L = np.array(
         [
-            [0.5, -0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0],
-            [-0.5, 0.5, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0],
-            [-0.5, -0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0],
-            [0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0],
-            [0.0, 0.0, 0.0, 0.0, 0.5, -0.35355339, 0.0, 0.0],
-            [0.0, 0.0, 0.0, 0.0, -0.35355339, 0.5, -0.35355339, -0.35355339],
-            [0.0, 0.0, 0.0, 0.0, 0.0, -0.35355339, 0.5, -0.5],
-            [0.0, 0.0, 0.0, 0.0, 0.0, -0.35355339, -0.5, 0.5],
+            [0.666667, -0.333333, -0.333333, -0.0, -0.0, -0.0, -0.0, -0.0],
+            [-0.333333, 0.666667, -0.333333, -0.0, -0.0, -0.0, -0.0, -0.0],
+            [-0.333333, -0.333333, 0.666667, -0.0, -0.0, -0.0, -0.0, -0.0],
+            [-0.0, -0.0, -0.0, 0.0, -0.0, -0.0, -0.0, -0.0],
+            [-0.0, -0.0, -0.0, -0.0, 0.5, -0.353553, -0.0, -0.0],
+            [-0.0, -0.0, -0.0, -0.0, -0.353553, 0.583333, -0.235702, -0.235702],
+            [-0.0, -0.0, -0.0, -0.0, -0.0, -0.235702, 0.666667, -0.333333],
+            [-0.0, -0.0, -0.0, -0.0, -0.0, -0.235702, -0.333333, 0.666667],
         ]
     )
     assert np.allclose(true_L, L2)
 
+    el_mwe = [
+        {1, 2, 3},
+        {1, 4, 5},
+        {1, 6, 7, 8},
+        {4, 9, 10, 11, 12},
+        {1, 13, 14, 15, 16},
+        {4, 17, 18},
+    ]
+    H_mwe = xgi.Hypergraph(el_mwe)
+    L_mwe = xgi.normalized_hypergraph_laplacian(H_mwe, sparse=False)
+    evals_mwe = eigvalsh(L_mwe)
+    assert np.all(evals_mwe >= 0)
+
+    # Weights error handling
+    ## Default
+    L_mwe_wtd = xgi.normalized_hypergraph_laplacian(H_mwe, weighted=True, sparse=False)
+    assert np.allclose(L_mwe, L_mwe_wtd)
+
+    ## Uniform weight
+    H_mwe.set_edge_attributes(2, name="weight")
+    L_mwe_wtd_uni = xgi.normalized_hypergraph_laplacian(
+        H_mwe, weighted=True, sparse=False
+    )
+    assert np.allclose(2 * L_mwe_wtd - np.eye(H_mwe.num_nodes), L_mwe_wtd_uni)
+
 
 def test_empty_order(edgelist6):
     H = xgi.Hypergraph(edgelist6)

xgi/linalg/laplacian_matrix.py

@@ -47,10 +47,15 @@
 from warnings import warn
 
 import numpy as np
-from scipy.sparse import csr_array, diags
+from scipy.sparse import csr_array, diags, identity
 
 from ..exception import XGIError
-from .hypergraph_matrix import adjacency_matrix, clique_motif_matrix, degree_matrix
+from .hypergraph_matrix import (
+    adjacency_matrix,
+    clique_motif_matrix,
+    degree_matrix,
+    incidence_matrix,
+)
 
 __all__ = [
     "laplacian",
@@ -183,13 +188,15 @@ def multiorder_laplacian(
     return (L_multi, rowdict) if index else L_multi
 
 
-def normalized_hypergraph_laplacian(H, sparse=True, index=False):
+def normalized_hypergraph_laplacian(H, weighted=False, sparse=True, index=False):
     """Compute the normalized Laplacian.
 
     Parameters
     ----------
     H : Hypergraph
         Hypergraph
+    weighted : bool, optional
+        whether or not to use hyperedge weightr, by default False (every edge weighted as 1).
     sparse : bool, optional
         whether or not the laplacian is sparse, by default True
     index : bool, optional
@@ -221,16 +228,39 @@ def normalized_hypergraph_laplacian(H, sparse=True, index=False):
             "Every node must be a member of an edge to avoid divide by zero error!"
         )
 
-    D = degree_matrix(H)
-    A, rowdict = clique_motif_matrix(H, sparse=sparse, index=True)
+    (
+        incidence,
+        rowdict,
+        _,  # Discard edge name-index mapping
+    ) = incidence_matrix(H, sparse=sparse, index=True)
+
+    Dv = degree_matrix(H)
+    De = np.sum(incidence, axis=0)
 
     if sparse:
-        Dinvsqrt = csr_array(diags(np.power(D, -0.5)))
+        Dv_invsqrt = csr_array(diags(np.power(Dv, -0.5)))
+        De_inv = diags(1 / De)
+
+        if weighted:
+            W = diags([H.edges[edge_idx].get("weight", 1) for edge_idx in H.edges])
+        else:
+            W = identity(H.num_edges)
+
         eye = csr_array((H.num_nodes, H.num_nodes))
         eye.setdiag(1)
     else:
-        Dinvsqrt = np.diag(np.power(D, -0.5))
+        Dv_invsqrt = np.diag(np.power(Dv, -0.5))
+        De_inv = np.diag(1 / De)
+
+        if weighted:
+            W = np.diag([H.edges[edge_idx].get("weight", 1) for edge_idx in H.edges])
+        else:
+            W = np.identity(H.num_edges)
+
         eye = np.eye(H.num_nodes)
 
-    L = 0.5 * (eye - Dinvsqrt @ A @ Dinvsqrt)
+    L = eye - (
+        Dv_invsqrt @ incidence @ W @ De_inv @ np.transpose(incidence) @ Dv_invsqrt
+    )
+
     return (L, rowdict) if index else L