"""Laplacian matrix of graphs.
All calculations here are done using the out-degree. For Laplacians using
in-degree, use `G.reverse(copy=False)` instead of `G` and take the transpose.
The `laplacian_matrix` function provides an unnormalized matrix,
while `normalized_laplacian_matrix`, `directed_laplacian_matrix`,
and `directed_combinatorial_laplacian_matrix` are all normalized.
"""
import networkx as nx
from networkx.utils import not_implemented_for
__all__ = [
"laplacian_matrix",
"normalized_laplacian_matrix",
"directed_laplacian_matrix",
"directed_combinatorial_laplacian_matrix",
"magnetic_laplacian_matrix",
]
[docs]
@nx._dispatchable(edge_attrs="weight")
def laplacian_matrix(G, nodelist=None, weight="weight"):
"""Returns the Laplacian matrix of G.
The graph Laplacian is the matrix L = D - A, where
A is the adjacency matrix and D is the diagonal matrix of node degrees.
Parameters
----------
G : graph
A NetworkX graph
nodelist : list, optional
The rows and columns are ordered according to the nodes in nodelist.
If nodelist is None, then the ordering is produced by G.nodes().
weight : string or None, optional (default='weight')
The edge data key used to compute each value in the matrix.
If None, then each edge has weight 1.
Returns
-------
L : SciPy sparse array
The Laplacian matrix of G.
Notes
-----
For MultiGraph, the edges weights are summed.
This returns an unnormalized matrix. For a normalized output,
use `normalized_laplacian_matrix`, `directed_laplacian_matrix`,
or `directed_combinatorial_laplacian_matrix`.
This calculation uses the out-degree of the graph `G`. To use the
in-degree for calculations instead, use `G.reverse(copy=False)` and
take the transpose.
See Also
--------
:func:`~networkx.convert_matrix.to_numpy_array`
normalized_laplacian_matrix
directed_laplacian_matrix
directed_combinatorial_laplacian_matrix
:func:`~networkx.linalg.spectrum.laplacian_spectrum`
Examples
--------
For graphs with multiple connected components, L is permutation-similar
to a block diagonal matrix where each block is the respective Laplacian
matrix for each component.
>>> G = nx.Graph([(1, 2), (2, 3), (4, 5)])
>>> print(nx.laplacian_matrix(G).toarray())
[[ 1 -1 0 0 0]
[-1 2 -1 0 0]
[ 0 -1 1 0 0]
[ 0 0 0 1 -1]
[ 0 0 0 -1 1]]
>>> edges = [
... (1, 2),
... (2, 1),
... (2, 4),
... (4, 3),
... (3, 4),
... ]
>>> DiG = nx.DiGraph(edges)
>>> print(nx.laplacian_matrix(DiG).toarray())
[[ 1 -1 0 0]
[-1 2 -1 0]
[ 0 0 1 -1]
[ 0 0 -1 1]]
Notice that node 4 is represented by the third column and row. This is because
by default the row/column order is the order of `G.nodes` (i.e. the node added
order -- in the edgelist, 4 first appears in (2, 4), before node 3 in edge (4, 3).)
To control the node order of the matrix, use the `nodelist` argument.
>>> print(nx.laplacian_matrix(DiG, nodelist=[1, 2, 3, 4]).toarray())
[[ 1 -1 0 0]
[-1 2 0 -1]
[ 0 0 1 -1]
[ 0 0 -1 1]]
This calculation uses the out-degree of the graph `G`. To use the
in-degree for calculations instead, use `G.reverse(copy=False)` and
take the transpose.
>>> print(nx.laplacian_matrix(DiG.reverse(copy=False)).toarray().T)
[[ 1 -1 0 0]
[-1 1 -1 0]
[ 0 0 2 -1]
[ 0 0 -1 1]]
References
----------
.. [1] Langville, Amy N., and Carl D. Meyer. Google’s PageRank and Beyond:
The Science of Search Engine Rankings. Princeton University Press, 2006.
"""
import scipy as sp
if nodelist is None:
nodelist = list(G)
A = nx.to_scipy_sparse_array(G, nodelist=nodelist, weight=weight, format="csr")
n, m = A.shape
D = sp.sparse.dia_array((A.sum(axis=1), 0), shape=(m, n)).tocsr()
return D - A
[docs]
@nx._dispatchable(edge_attrs="weight")
def normalized_laplacian_matrix(G, nodelist=None, weight="weight"):
r"""Returns the normalized Laplacian matrix of G.
The normalized graph Laplacian is the matrix
.. math::
N = D^{-1/2} L D^{-1/2}
where `L` is the graph Laplacian and `D` is the diagonal matrix of
node degrees [1]_.
Parameters
----------
G : graph
A NetworkX graph
nodelist : list, optional
The rows and columns are ordered according to the nodes in nodelist.
If nodelist is None, then the ordering is produced by G.nodes().
weight : string or None, optional (default='weight')
The edge data key used to compute each value in the matrix.
If None, then each edge has weight 1.
Returns
-------
N : SciPy sparse array
The normalized Laplacian matrix of G.
Notes
-----
For MultiGraph, the edges weights are summed.
See :func:`to_numpy_array` for other options.
If the Graph contains selfloops, D is defined as ``diag(sum(A, 1))``, where A is
the adjacency matrix [2]_.
This calculation uses the out-degree of the graph `G`. To use the
in-degree for calculations instead, use `G.reverse(copy=False)` and
take the transpose.
For an unnormalized output, use `laplacian_matrix`.
Examples
--------
>>> import numpy as np
>>> edges = [
... (1, 2),
... (2, 1),
... (2, 4),
... (4, 3),
... (3, 4),
... ]
>>> DiG = nx.DiGraph(edges)
>>> print(nx.normalized_laplacian_matrix(DiG).toarray())
[[ 1. -0.70710678 0. 0. ]
[-0.70710678 1. -0.70710678 0. ]
[ 0. 0. 1. -1. ]
[ 0. 0. -1. 1. ]]
Notice that node 4 is represented by the third column and row. This is because
by default the row/column order is the order of `G.nodes` (i.e. the node added
order -- in the edgelist, 4 first appears in (2, 4), before node 3 in edge (4, 3).)
To control the node order of the matrix, use the `nodelist` argument.
>>> print(nx.normalized_laplacian_matrix(DiG, nodelist=[1, 2, 3, 4]).toarray())
[[ 1. -0.70710678 0. 0. ]
[-0.70710678 1. 0. -0.70710678]
[ 0. 0. 1. -1. ]
[ 0. 0. -1. 1. ]]
>>> G = nx.Graph(edges)
>>> print(nx.normalized_laplacian_matrix(G).toarray())
[[ 1. -0.70710678 0. 0. ]
[-0.70710678 1. -0.5 0. ]
[ 0. -0.5 1. -0.70710678]
[ 0. 0. -0.70710678 1. ]]
See Also
--------
laplacian_matrix
normalized_laplacian_spectrum
directed_laplacian_matrix
directed_combinatorial_laplacian_matrix
References
----------
.. [1] Fan Chung-Graham, Spectral Graph Theory,
CBMS Regional Conference Series in Mathematics, Number 92, 1997.
.. [2] Steve Butler, Interlacing For Weighted Graphs Using The Normalized
Laplacian, Electronic Journal of Linear Algebra, Volume 16, pp. 90-98,
March 2007.
.. [3] Langville, Amy N., and Carl D. Meyer. Google’s PageRank and Beyond:
The Science of Search Engine Rankings. Princeton University Press, 2006.
"""
import numpy as np
import scipy as sp
if nodelist is None:
nodelist = list(G)
A = nx.to_scipy_sparse_array(G, nodelist=nodelist, weight=weight, format="csr")
n, _ = A.shape
diags = A.sum(axis=1)
D = sp.sparse.dia_array((diags, 0), shape=(n, n)).tocsr()
L = D - A
with np.errstate(divide="ignore"):
diags_sqrt = 1.0 / np.sqrt(diags)
diags_sqrt[np.isinf(diags_sqrt)] = 0
DH = sp.sparse.dia_array((diags_sqrt, 0), shape=(n, n)).tocsr()
return DH @ (L @ DH)
[docs]
@not_implemented_for("multigraph")
@not_implemented_for("undirected")
@nx._dispatchable(edge_attrs="weight")
def magnetic_laplacian_matrix(
G, *, nodelist=None, normalized=False, q=0.25, weight="weight"
):
r"""Returns the magnetic Laplacian matrix of DiGraph G.
The magnetic Laplacian matrix (also called the q-magnetic Laplacian) is a
Hermitian matrix for directed graphs that encodes edge directionality
via complex phases [1]_.
Given a weighted directed graph, :math:`G = (V, E, W)`, with :math:`W`
the weighted adjacency matrix, the symmetrized weighted adjacency
matrix is defined as :math:`W' = 0.5 (W + W^{T})`. A skew-symmetric term
:math:`\delta` is introduced to encode directionality, where
..math::
\delta_{jk} = \begin{cases}
+1 & \text{if} j \to k \text{ is an edge and } k \to j \text{ is not}, \\
-1 & \text{if} k \to j \text{ is an edge and } j \to k \text{ is not}, \\
0 & \text{if both or neither edge is present.}
\end{cases}
Then, the magnetic Laplacian matrix is defined as:
.. math::
L^{(q)} := D - H^{(q)}
where :math:`H^{(q)}` is the Hermitian adjacency matrix with entries
:math:`H^{(q)}_{jk} = W'_{jk} e^{2\pi i q \delta_{jk}}`, and :math:`D` is the
degree matrix associated with the symmetrized weight adjacency matrix :math:`W'`.
If `normalized` is True, compute the normalized version using the Moore-Penrose
inverse :math:`D^{+}}` of the degree matrix :math:`D`. The normalized formula
is then:
.. math::
L^{(q)}_{norm} = (D^{+})^{0.5} L^{(q)} (D^{+})^{0.5}
Parameters
----------
G : DiGraph
A directed graph
nodelist : list, optional (default=list(G))
Node ordering for row/columns.
normalized : bool, optional (default=False)
Bool that encodes if return the magnetic Laplacian or the normalized magnetic
Laplacian. If True returns the normalized version.
q : float, optional (default=0.25)
The phase of the magnetic potential is the charge parameter 0 <= q <= 0.5.
At q=0 returns the standard Laplacian.
weight : string or None, optional (default='weight')
Edge attribute key for weights. If None, all edges have weight 1.
Returns
-------
L : SciPy sparse array (complex dtype)
The magnetic Laplacian matrix of `G` if not `normalized`
and the normalized version if `normalized`
Raises
------
ValueError
If q is not between 0 and 0.5
NetworkXNotImplemented
If `G` is undirected or a multigraph
References
----------
.. [1] Fanuel, M., Alaíz, C. M., Fernández, Á., & Suykens, J. A. (2018).
Magnetic eigenmaps for the visualization of directed graphs.
Applied and Computational Harmonic Analysis, 44(1), 189–199.
<https://doi.org/10.1016/j.acha.2017.01.004>
"""
from collections import defaultdict
import numpy as np
import scipy as sp
if nodelist is None:
nodelist = list(G)
# Build Hermitian adjacency H
n = len(nodelist)
node_index = {v: i for i, v in enumerate(nodelist)}
if not (0 <= q <= 0.5):
raise ValueError("Parameter q must be a value between 0 and 0.5")
phase = 2 * np.pi * q
# Find "phase" matrix
# Dict to encode where have been added symetries
delta_phase_edge = defaultdict(int)
matrix_weights = defaultdict(float)
phases = {1: np.exp(1j * phase), -1: np.exp(-1j * phase), 0: 1}
for u, v, wt in G.edges(data=weight, default=1):
if u not in node_index or v not in node_index:
continue
ui, vi = node_index[u], node_index[v]
delta_phase_edge[(ui, vi)] += 1
delta_phase_edge[(vi, ui)] -= 1
matrix_weights[(ui, vi)] += 0.5 * wt
matrix_weights[(vi, ui)] += 0.5 * wt
rows, cols, data = [], [], []
for u, v in G.edges():
if u not in node_index or v not in node_index:
continue
ui, vi = node_index[u], node_index[v]
if ui != vi:
if delta_phase_edge[(ui, vi)] == 0:
rows.append(ui)
cols.append(vi)
data.append(matrix_weights[(ui, vi)])
else:
rows.append(ui)
cols.append(vi)
data.append(
matrix_weights[(ui, vi)] * phases[delta_phase_edge[(ui, vi)]]
)
rows.append(vi)
cols.append(ui)
data.append(
matrix_weights[(vi, ui)] * phases[delta_phase_edge[(vi, ui)]]
)
H = sp.sparse.csr_array((data, (rows, cols)), shape=(n, n), dtype=complex)
# Build degree matrix D
diags = np.abs(H).sum(axis=1).ravel()
if normalized:
with np.errstate(divide="ignore"):
diags_sqrt = 1.0 / np.sqrt(diags)
diags_sqrt[np.isinf(diags_sqrt)] = 0
DH = sp.sparse.dia_array((diags_sqrt, 0), shape=(n, n)).tocsr()
H = DH @ (H @ DH)
diags = np.ones(
n,
)
D = sp.sparse.dia_array((diags, 0), shape=(n, n), dtype=complex).tocsr()
return D - H
###############################################################################
# Code based on work from https://github.com/bjedwards
[docs]
@not_implemented_for("undirected")
@not_implemented_for("multigraph")
@nx._dispatchable(edge_attrs="weight")
def directed_laplacian_matrix(
G, nodelist=None, weight="weight", walk_type=None, alpha=0.95
):
r"""Returns the directed Laplacian matrix of G.
The graph directed Laplacian is the matrix
.. math::
L = I - \frac{1}{2} \left (\Phi^{1/2} P \Phi^{-1/2} + \Phi^{-1/2} P^T \Phi^{1/2} \right )
where `I` is the identity matrix, `P` is the transition matrix of the
graph, and `\Phi` a matrix with the Perron vector of `P` in the diagonal and
zeros elsewhere [1]_.
Depending on the value of walk_type, `P` can be the transition matrix
induced by a random walk, a lazy random walk, or a random walk with
teleportation (PageRank).
Parameters
----------
G : DiGraph
A NetworkX graph
nodelist : list, optional
The rows and columns are ordered according to the nodes in nodelist.
If nodelist is None, then the ordering is produced by G.nodes().
weight : string or None, optional (default='weight')
The edge data key used to compute each value in the matrix.
If None, then each edge has weight 1.
walk_type : string or None, optional (default=None)
One of ``"random"``, ``"lazy"``, or ``"pagerank"``. If ``walk_type=None``
(the default), then a value is selected according to the properties of `G`:
- ``walk_type="random"`` if `G` is strongly connected and aperiodic
- ``walk_type="lazy"`` if `G` is strongly connected but not aperiodic
- ``walk_type="pagerank"`` for all other cases.
alpha : real
(1 - alpha) is the teleportation probability used with pagerank
Returns
-------
L : NumPy matrix
Normalized Laplacian of G.
Notes
-----
Only implemented for DiGraphs
The result is always a symmetric matrix.
This calculation uses the out-degree of the graph `G`. To use the
in-degree for calculations instead, use `G.reverse(copy=False)` and
take the transpose.
See Also
--------
laplacian_matrix
normalized_laplacian_matrix
directed_combinatorial_laplacian_matrix
References
----------
.. [1] Fan Chung (2005).
Laplacians and the Cheeger inequality for directed graphs.
Annals of Combinatorics, 9(1), 2005
"""
import numpy as np
import scipy as sp
# NOTE: P has type ndarray if walk_type=="pagerank", else csr_array
P = _transition_matrix(
G, nodelist=nodelist, weight=weight, walk_type=walk_type, alpha=alpha
)
n, m = P.shape
evals, evecs = sp.sparse.linalg.eigs(P.T, k=1)
v = evecs.flatten().real
p = v / v.sum()
# p>=0 by Perron-Frobenius Thm. Use abs() to fix roundoff across zero gh-6865
sqrtp = np.sqrt(np.abs(p))
Q = (
sp.sparse.dia_array((sqrtp, 0), shape=(n, n)).tocsr()
[docs]
@ P
@ sp.sparse.dia_array((1.0 / sqrtp, 0), shape=(n, n)).tocsr()
)
# NOTE: This could be sparsified for the non-pagerank cases
I = np.identity(len(G))
return I - (Q + Q.T) / 2.0
@not_implemented_for("undirected")
@not_implemented_for("multigraph")
@nx._dispatchable(edge_attrs="weight")
def directed_combinatorial_laplacian_matrix(
G, nodelist=None, weight="weight", walk_type=None, alpha=0.95
):
r"""Return the directed combinatorial Laplacian matrix of G.
The graph directed combinatorial Laplacian is the matrix
.. math::
L = \Phi - \frac{1}{2} \left (\Phi P + P^T \Phi \right)
where `P` is the transition matrix of the graph and `\Phi` a matrix
with the Perron vector of `P` in the diagonal and zeros elsewhere [1]_.
Depending on the value of walk_type, `P` can be the transition matrix
induced by a random walk, a lazy random walk, or a random walk with
teleportation (PageRank).
Parameters
----------
G : DiGraph
A NetworkX graph
nodelist : list, optional
The rows and columns are ordered according to the nodes in nodelist.
If nodelist is None, then the ordering is produced by G.nodes().
weight : string or None, optional (default='weight')
The edge data key used to compute each value in the matrix.
If None, then each edge has weight 1.
walk_type : string or None, optional (default=None)
One of ``"random"``, ``"lazy"``, or ``"pagerank"``. If ``walk_type=None``
(the default), then a value is selected according to the properties of `G`:
- ``walk_type="random"`` if `G` is strongly connected and aperiodic
- ``walk_type="lazy"`` if `G` is strongly connected but not aperiodic
- ``walk_type="pagerank"`` for all other cases.
alpha : real
(1 - alpha) is the teleportation probability used with pagerank
Returns
-------
L : NumPy matrix
Combinatorial Laplacian of G.
Notes
-----
Only implemented for DiGraphs
The result is always a symmetric matrix.
This calculation uses the out-degree of the graph `G`. To use the
in-degree for calculations instead, use `G.reverse(copy=False)` and
take the transpose.
See Also
--------
laplacian_matrix
normalized_laplacian_matrix
directed_laplacian_matrix
References
----------
.. [1] Fan Chung (2005).
Laplacians and the Cheeger inequality for directed graphs.
Annals of Combinatorics, 9(1), 2005
"""
import scipy as sp
P = _transition_matrix(
G, nodelist=nodelist, weight=weight, walk_type=walk_type, alpha=alpha
)
n, m = P.shape
evals, evecs = sp.sparse.linalg.eigs(P.T, k=1)
v = evecs.flatten().real
p = v / v.sum()
# NOTE: could be improved by not densifying
Phi = sp.sparse.dia_array((p, 0), shape=(n, n)).toarray()
return Phi - (Phi @ P + P.T @ Phi) / 2.0
def _transition_matrix(G, nodelist=None, weight="weight", walk_type=None, alpha=0.95):
"""Returns the transition matrix of G.
This is a row stochastic giving the transition probabilities while
performing a random walk on the graph. Depending on the value of walk_type,
P can be the transition matrix induced by a random walk, a lazy random walk,
or a random walk with teleportation (PageRank).
Parameters
----------
G : DiGraph
A NetworkX graph
nodelist : list, optional
The rows and columns are ordered according to the nodes in nodelist.
If nodelist is None, then the ordering is produced by G.nodes().
weight : string or None, optional (default='weight')
The edge data key used to compute each value in the matrix.
If None, then each edge has weight 1.
walk_type : string or None, optional (default=None)
One of ``"random"``, ``"lazy"``, or ``"pagerank"``. If ``walk_type=None``
(the default), then a value is selected according to the properties of `G`:
- ``walk_type="random"`` if `G` is strongly connected and aperiodic
- ``walk_type="lazy"`` if `G` is strongly connected but not aperiodic
- ``walk_type="pagerank"`` for all other cases.
alpha : real
(1 - alpha) is the teleportation probability used with pagerank
Returns
-------
P : numpy.ndarray
transition matrix of G.
Raises
------
NetworkXError
If walk_type not specified or alpha not in valid range
"""
import numpy as np
import scipy as sp
if walk_type is None:
if nx.is_strongly_connected(G):
if nx.is_aperiodic(G):
walk_type = "random"
else:
walk_type = "lazy"
else:
walk_type = "pagerank"
A = nx.to_scipy_sparse_array(G, nodelist=nodelist, weight=weight, dtype=float)
n, m = A.shape
if walk_type in ["random", "lazy"]:
DI = sp.sparse.dia_array((1.0 / A.sum(axis=1), 0), shape=(n, n)).tocsr()
if walk_type == "random":
P = DI @ A
else:
I = sp.sparse.eye_array(n, format="csr")
P = (I + DI @ A) / 2.0
elif walk_type == "pagerank":
if not (0 < alpha < 1):
raise nx.NetworkXError("alpha must be between 0 and 1")
# this is using a dense representation. NOTE: This should be sparsified!
A = A.toarray()
# add constant to dangling nodes' row
A[A.sum(axis=1) == 0, :] = 1 / n
# normalize
A = A / A.sum(axis=1)[np.newaxis, :].T
P = alpha * A + (1 - alpha) / n
else:
raise nx.NetworkXError("walk_type must be random, lazy, or pagerank")
return P
