diff --git a/tests/linalg/test_matrix.py b/tests/linalg/test_matrix.py index cd59ff2f6..07e44379e 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): @@ -605,24 +607,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): 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