Repository: xgi-org/xgi

Issue:
Bug: Negative eigenvalues from `normalized_hypergraph_laplacian`
# Summary
The normalized hypergraph Laplacian of a hypergraph with non-uniformly sized edges is capable of returning negative eigenvalues.

# Example
Below is a (minimal?) working example
```python
import xgi
import numpy as np
from scipy.linalg import eigh

E = [
    {1, 2, 3},
    {1, 4, 5},
    {1, 6, 7, 8},
    {4, 9, 10, 11, 12},
    {1, 13, 14, 15, 16},
    {4, 17, 18}
]
H = xgi.Hypergraph(E)

L = xgi.normalized_hypergraph_laplacian(H, sparse=False)
evals, evecs = eigh(L)

assert np.all(evals >= 0)
```

# Suspected problem 
The [referenced paper](https://proceedings.neurips.cc/paper/2006/hash/dff8e9c2ac33381546d96deea9922999-Abstract.html) [1] in the implementation states that the Laplacian should have a trivial eigenvalue, eigenvector pair with eigenvalue 0 and eigenvector $(\sqrt{k_i})_i$. Since the Laplacian is real symmetric and positive semi-definite, it should also have non-negative real eigenvalues so I believe that eigenvalues are impossible.

Looking through the [source implementation](https://github.com/xgi-org/xgi/blob/main/xgi/linalg/laplacian_matrix.py#L186), I suspect the problem is that the current implementation utilizes a simplification of the Laplacian that is only valid in 2-uniform hypergraphs as derived in [1].

# Suggested fix
I have been able to circumvent this error by relying on the definition directly in [1], i.e., 
$$\mathcal{L} = I - D_v^{-\frac{1}{2}}HD_e^{-1}H^{T}D_v^{-\frac{1}{2}}$$ (for unweighted hypergraphs).

I am happy to create a pull request!


Also, this is my first time submitting an issue to a public repo - sorry if I have done anything wrong with it!

---

[1] “Learning with Hypergraphs: Clustering, Classification, and Embedding” by Dengyong Zhou, Jiayuan Huang, Bernhard Schölkopf Advances in Neural Information Processing Systems (2006)

Discussion:
@maximelucas: Hi @kaiser-dan, thanks for raising this issue, you exposed things very clearly!

After a quick check of the ref, I agree with you, it looks like we should've used the general formula. 

For future reference, the original ref even works for weighted hypergraphs by adding a weight matrix $W$, where $\Delta$ is the Laplacian (top of page 5):

![Image](https://github.com/user-attachments/assets/412cbdf1-8e2a-4494-af42-a51508d8d40e)

Feel free to create a pull request with your current solution, thanks so much!
@kaiser-dan: Okie done, I have created a pull request #648

Solution patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 610dddc64..564f7e527 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -13,7 +13,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
+        python-version: ['3.10', '3.11', '3.12', '3.13']
     steps:
     - uses: actions/checkout@v4
 
@@ -112,7 +112,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu, macos, windows]
-        python-version: ['3.9', '3.10', '3.11', '3.12', '3.13']
+        python-version: ['3.10', '3.11', '3.12', '3.13']
     steps:
     - uses: actions/checkout@v4
 
diff --git a/README.md b/README.md
index 699d81374..8d8565165 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 <img src="https://raw.githubusercontent.com/xgi-org/xgi/main/logo/logo.svg" alt="XGI" width="50%"></img>
 
 [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
-[![Supports Python versions 3.9 and above.](https://img.shields.io/badge/Python%20versions%20supported-3.9%2B-forest
+[![Supports Python versions 3.10 and above.](https://img.shields.io/badge/Python%20versions%20supported-3.10%2B-forest
 )](https://www.repostatus.org/#active)
 [![Test Status](https://github.com/xgi-org/xgi/workflows/test/badge.svg?branch=main)](https://github.com/xgi-org/xgi/actions?query=workflow%3A%22test%22)
 [![codecov](https://codecov.io/gh/xgi-org/xgi/branch/main/graph/badge.svg?token=BI6TX2WDSG)](https://codecov.io/gh/xgi-org/xgi)
@@ -43,7 +43,7 @@ Comple**X** **G**roup **I**nteractions (**XGI**) is a Python package for higher-
 * Draw higher-order networks in a variety of **visually striking ways** (See our [gallery](https://xgi.readthedocs.io/en/stable/auto_examples/index.html) for several examples.)
 
 ## Installation<a id="installation"></a>
-XGI runs on Python 3.9 or higher.
+XGI runs on Python 3.10 or higher.
 
 To install the latest version of XGI, run the following command:
 ```sh
diff --git a/pyproject.toml b/pyproject.toml
index d39a01a27..e5e6cc12d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,7 @@ name = "xgi"
 dynamic = ["version"]
 description = "XGI is a Python package for higher-order networks."
 readme = "README.md"
-requires-python = ">=3.9"
+requires-python = ">=3.10"
 license = {file = "LICENSE.md"}
 authors = [
     {name = "XGI Developers", email = "nicholas.landry@virginia.edu"}
@@ -30,7 +30,6 @@ classifiers = [
     "License :: OSI Approved :: BSD License",
     "Operating System :: OS Independent",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
@@ -43,8 +42,8 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Physics",
 ]
 dependencies = [
-    "numpy>=1.20",
-    "scipy>=1.8",
+    "numpy>=1.23.5",
+    "scipy>=1.15",
     "pandas>=1.3",
     "networkx>=2.7",
     "requests>=2.0",
diff --git a/requirements/default.txt b/requirements/default.txt
index c25b72295..0b67041e9 100644
--- a/requirements/default.txt
+++ b/requirements/default.txt
@@ -1,7 +1,7 @@
 # Generated via tools/generate_requirements.py and pre-commit hook.
 # Do not edit this file; modify pyproject.toml instead.
-numpy>=1.20
-scipy>=1.8
+numpy>=1.23.5
+scipy>=1.15
 pandas>=1.3
 networkx>=2.7
 requests>=2.0
diff --git a/tools/generate_requirements.py b/tools/generate_requirements.py
index 35fbe8a16..304defc9f 100644
--- a/tools/generate_requirements.py
+++ b/tools/generate_requirements.py
@@ -1,6 +1,8 @@
 #!/usr/bin/env python
 """Generate requirements/*.txt files from pyproject.toml.
 
+Run from repo root folder.
+
 Copied from NetworkX: https://github.com/networkx/networkx/blob/main/tools/generate_requirements.py
 """
 
diff --git a/xgi/linalg/laplacian_matrix.py b/xgi/linalg/laplacian_matrix.py
index 1ebdd686a..47a3ff323 100644
--- a/xgi/linalg/laplacian_matrix.py
+++ b/xgi/linalg/laplacian_matrix.py
@@ -47,10 +47,14 @@
 from warnings import warn
 
 import numpy as np
-from scipy.sparse import csr_array, diags
+from scipy.sparse import csr_array, diags_array, eye_array
 
 from ..exception import XGIError
-from .hypergraph_matrix import adjacency_matrix, clique_motif_matrix, degree_matrix
+from .hypergraph_matrix import (
+    adjacency_matrix,
+    degree_matrix,
+    incidence_matrix,
+)
 
 __all__ = [
     "laplacian",
@@ -102,7 +106,7 @@ def laplacian(H, order=1, sparse=False, rescale_per_node=False, index=False):
         return (L, {}) if index else L
 
     if sparse:
-        K = csr_array(diags(degree_matrix(H, order=order)))
+        K = diags_array(degree_matrix(H, order=order), format="csr")
     else:
         K = np.diag(degree_matrix(H, order=order))
 
@@ -183,13 +187,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 weights, 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 +227,31 @@ 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 weighted:
+        weights = [H.edges[edge_idx].get("weight", 1) for edge_idx in H.edges]
+    else:
+        weights = [1] * H.num_edges
 
     if sparse:
-        Dinvsqrt = csr_array(diags(np.power(D, -0.5)))
-        eye = csr_array((H.num_nodes, H.num_nodes))
-        eye.setdiag(1)
+        Dv_invsqrt = diags_array(np.power(Dv, -0.5), format="csr")
+        De_inv = diags_array(1 / De, format="csr")
+        W = diags_array(weights, format="csr")
+        eye = eye_array(H.num_nodes, format="csr")
     else:
-        Dinvsqrt = np.diag(np.power(D, -0.5))
+        Dv_invsqrt = np.diag(np.power(Dv, -0.5))
+        De_inv = np.diag(1 / De)
+        W = np.diag(weights)
         eye = np.eye(H.num_nodes)
 
-    L = 0.5 * (eye - Dinvsqrt @ A @ Dinvsqrt)
+    L = eye - Dv_invsqrt @ incidence @ W @ De_inv @ incidence.T @ Dv_invsqrt
+
     return (L, rowdict) if index else L


Test patch:
diff --git a/tests/linalg/test_matrix.py b/tests/linalg/test_matrix.py
index cd59ff2f6..75e6f8b7c 100644
--- a/tests/linalg/test_matrix.py
+++ b/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):
@@ -255,6 +257,23 @@ def test_adjacency_matrix(edgelist1, edgelist4):
     assert np.all(A7 == A7_sp.todense())
 
 
+def test_fix_649(): # make sure diagonal is fully zero when sparse
+
+    H = xgi.load_xgi_data(dataset="email-enron")
+    H.cleanup(isolates=True, singletons=False, connected=True, relabel=False, multiedges=False)
+
+    weighted = True
+    order = 1
+
+    adj = xgi.adjacency_matrix(H, order=order, sparse=False, weighted=weighted, index=False)
+
+    adj_sparse = xgi.adjacency_matrix(H, order=order, sparse=True, weighted=weighted, index=False)
+    adj_sparse = adj_sparse.todense()
+
+    assert np.all(np.diag(adj_sparse)==0)
+
+
+
 def test_laplacian(edgelist2, edgelist6):
     el1 = edgelist6
     H1 = xgi.Hypergraph(el1)
@@ -605,24 +624,64 @@ 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)
+    # Eigenvalues are all non-negative
+    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}
-    true_L = np.array(
+    assert np.allclose(L3.toarray(), L2)
+
+    # 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)
+
+    # Exact Laplacian calculation
+    true_L3 = 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)
+    assert np.allclose(true_L3, L3.toarray())
+
+
+def test_fix_647():
+    el = [
+        {1, 2, 3},
+        {1, 4, 5},
+        {1, 6, 7, 8},
+        {4, 9, 10, 11, 12},
+        {1, 13, 14, 15, 16},
+        {4, 17, 18},
+    ]
+    H = xgi.Hypergraph(el)
+    L = xgi.normalized_hypergraph_laplacian(H, sparse=False)
+
+    # Eigenvalues non-negative
+    evals_mwe = eigvalsh(L)
+    assert np.all(evals_mwe >= 0)
+
+    # Weights error handling
+    ## Default value when "weight" attribute unavailable
+    L_wtd = xgi.normalized_hypergraph_laplacian(H, weighted=True, sparse=False)
+    assert np.allclose(L, L_wtd)
+
+    ## Uniform weight
+    H_wtd = H.copy()
+    H_wtd.set_edge_attributes(2, name="weight")
+    L_wtd_uniform = xgi.normalized_hypergraph_laplacian(
+        H_wtd, weighted=True, sparse=False
+    )
+    assert np.allclose(2 * L_wtd - np.eye(H_wtd.num_nodes), L_wtd_uniform)
 
 
 def test_empty_order(edgelist6):