Repository: networkx/networkx

Issue:
Unstable and Incorrect "nx.eigenvector_centrality_numpy" result when the graph is disconnected
Hello! Sorry for bothering you again!
I found that when the graph is disconnected, the result of the `nx.eigenvector_centrality_numpy` is unstable and different from the `nx.eigenvector_centrality`. Here is the reduced test case with 4 nodes and 2 edges:

### Step to Reproduce:
```Python
import networkx as nx

G = nx.from_edgelist([(1, 2), (3, 4)])
print(nx.eigenvector_centrality_numpy(G))
G = nx.from_edgelist([(1, 2), (3, 4)])
print(nx.eigenvector_centrality_numpy(G))
print(nx.eigenvector_centrality(G))
```
#### Results:
```
#nx.eigenvector_centrality_numpy (unstable!)
{1: -0.14904837055016898, 2: -0.14904837055016898, 3: 0.6912196345853752, 4: 0.6912196345853752}
{1: -0.18977204610658815, 2: -0.1897720461065881, 3: 0.6811655969854312, 4: 0.6811655969854313}
#nx.eigenvector_centrality
{1: 0.5, 2: 0.5, 3: 0.5, 4: 0.5}
```
I believe that it may be related to a logic bug in the implementation of `nx.eigenvector_centrality_numpy`. Would it be possible for you to further confirm and investigate it?

Best regards,
Joye

### Environments
NetworkX: 3.1
Python: 3.10


Discussion:
Peiffap: All sensible suggestions, thanks for taking a look!

It'd definitely be nice to have benchmarks for this kind of change (re: what I said in #7527). I could add a simple benchmark right now, but I think there might have to be a reorganisation of how benchmarking is done at some point in the near future if performance is something you want to prioritize.
Peiffap: I added a simple benchmarking function. I **didn't test this locally** (it's a pain getting it to work without conda or mamba installed...), so we should make sure it's working as expected before merging
eriknw: FYI: somewhat related, we have, and are, exploring caching results of functions (see #7283). For example, this can allow us to compute connectedness of a graph only once, which may help mitigate the performance penalty of checking connectedness (and other such properties) here and elsewhere.

At PyConUS @dschult, @rlratzel, @MridulS, and I discussed that caching of function results should be added _narrowly_, _incrementally_, and _as use cases arise_, unlike what was done in #7283. Hence maybe caching could begin with e.g. `is_connected`, `is_strongly_connected`, `is_weakly_connected`, `number_connected_components`, `number_strongly_connected_components`, `number_weakly_connected_components`, or just some of these.

This doesn't change anything that needs done in this PR.  It seems like a good idea to me to check for properties such as connectedness to ensure that an algorithm is well-defined for the input graph.
rossbar: I went ahead and split up the benchmark in a naive way in order to unstick this - maybe @MridulS can doublecheck that this makes sense. Anyways, my approval of the changes themselves (regardless of benchmarking changes) stands!

Patch:
diff --git a/benchmarks/benchmarks/benchmark_algorithms.py b/benchmarks/benchmarks/benchmark_algorithms.py
index 8991912c7..34df6ed00 100644
--- a/benchmarks/benchmarks/benchmark_algorithms.py
+++ b/benchmarks/benchmarks/benchmark_algorithms.py
@@ -48,3 +48,28 @@ class AlgorithmBenchmarks:
 
     def time_average_clustering(self, graph):
         _ = nx.average_clustering(self.graphs_dict[graph])
+
+
+class AlgorithmBenchmarksConnectedGraphsOnly:
+    timeout = 120
+    nodes = 100
+    _graphs = [
+        nx.erdos_renyi_graph(nodes, 0.1),
+        nx.erdos_renyi_graph(nodes, 0.5),
+        nx.erdos_renyi_graph(nodes, 0.9),
+    ]
+    params = [
+        "Erdos Renyi (100, 0.1)",
+        "Erdos Renyi (100, 0.5)",
+        "Erdos Renyi (100, 0.9)",
+    ]
+
+    param_names = ["graph"]
+
+    def setup(self, graph):
+        self.graphs_dict = dict(zip(self.params, self._graphs))
+
+    def time_eigenvector_centrality_numpy(self, graph):
+        # Added to ensure the connectivity check doesn't affect
+        # performance too much (see gh-6888, gh-7549).
+        _ = nx.eigenvector_centrality_numpy(self.graphs_dict[graph])
diff --git a/doc/developer/about_us.rst b/doc/developer/about_us.rst
index e9e41adf2..34fadeb5d 100644
--- a/doc/developer/about_us.rst
+++ b/doc/developer/about_us.rst
@@ -154,6 +154,7 @@ to add your name to the bottom of the list.
 - Aaron Zolnai-Lucas, GitHub: `aaronzo <https://github.com/aaronzo>`_, LinkedIn: `aaronzolnailucas <https://www.linkedin.com/in/aaronzolnailucas/>`_
 - Erik Welch, GitHub: `eriknw <https://github.com/eriknw>`_, LinkedIn: `eriknwelch <https://www.linkedin.com/in/eriknwelch/>`_
 - Mohamed Rezk, Github: `mohamedrezk122 <https://github.com/mohamedrezk122>`_
+- Orion Sehn, GitHub: `OrionSehn <https://github.com/OrionSehn>`_ 
 
 A supplementary (but still incomplete) list of contributors is given by the
 list of names that have commits in ``networkx``'s
diff --git a/doc/developer/projects.rst b/doc/developer/projects.rst
index fbd3b4e2c..b7ef0e714 100644
--- a/doc/developer/projects.rst
+++ b/doc/developer/projects.rst
@@ -112,6 +112,11 @@ Centrality Atlas
 Completed Projects
 ==================
 
+- `Revisiting and expanding nx-parallel`_
+    - Program: Google Summer of Code 2024
+    - Contributor: `@Schefflera-Arboricola <https://github.com/Schefflera-Arboricola>`__
+    - Link to Proposal: `GSoC 2024: Revisiting and expanding nx-parallel <https://github.com/networkx/archive/blob/main/proposals-gsoc/GSoC-2024-Revisiting-and-expanding-nx-parallel.pdf>`_
+
 - `VF2++ algorithm for graph isomorphism`_
     - Program: Google Summer of Code 2022
     - Contributor: `@kpetridis24 <https://github.com/kpetridis24/>`__
@@ -145,6 +150,7 @@ Completed Projects
     - Contributor: `@MridulS <https://github.com/MridulS/>`__
     - Link to Proposal: `GSoC 2015: NetworkX 2.0 API <https://github.com/networkx/archive/blob/main/proposals-gsoc/GSoC-2015-NetworkX-2.0-api.md>`__
 
+.. _`Revisiting and expanding nx-parallel`: https://github.com/Schefflera-Arboricola/blogs/tree/main/networkx/GSoC24
 .. _`VF2++ algorithm for graph isomorphism`: https://github.com/networkx/networkx/pull/5788
 .. _`Louvain community detection algorithm`: https://github.com/networkx/networkx/pull/4929
 .. _`Asadpour algorithm for directed travelling salesman problem`: https://github.com/networkx/networkx/pull/4740
diff --git a/doc/reference/algorithms/distance_measures.rst b/doc/reference/algorithms/distance_measures.rst
index 172e5b6dd..00da08c10 100644
--- a/doc/reference/algorithms/distance_measures.rst
+++ b/doc/reference/algorithms/distance_measures.rst
@@ -9,6 +9,7 @@ Distance Measures
    barycenter
    center
    diameter
+   harmonic_diameter
    eccentricity
    effective_graph_resistance
    kemeny_constant
diff --git a/networkx/algorithms/centrality/eigenvector.py b/networkx/algorithms/centrality/eigenvector.py
index 2cef241b0..b8cf63e8d 100644
--- a/networkx/algorithms/centrality/eigenvector.py
+++ b/networkx/algorithms/centrality/eigenvector.py
@@ -196,7 +196,7 @@ def eigenvector_centrality(G, max_iter=100, tol=1.0e-6, nstart=None, weight=None
 
 @nx._dispatchable(edge_attrs="weight")
 def eigenvector_centrality_numpy(G, weight=None, max_iter=50, tol=0):
-    r"""Compute the eigenvector centrality for the graph G.
+    r"""Compute the eigenvector centrality for the graph `G`.
 
     Eigenvector centrality computes the centrality for a node by adding
     the centrality of its predecessors. The centrality for node $i$ is the
@@ -208,7 +208,7 @@ def eigenvector_centrality_numpy(G, weight=None, max_iter=50, tol=0):
 
          \lambda x^T = x^T A,
 
-    where $A$ is the adjacency matrix of the graph G. By definition of
+    where $A$ is the adjacency matrix of the graph `G`. By definition of
     row-column product, the equation above is equivalent to
 
     .. math::
@@ -220,37 +220,42 @@ def eigenvector_centrality_numpy(G, weight=None, max_iter=50, tol=0):
     $\lambda$. In the case of undirected graphs, $x$ also solves the familiar
     right-eigenvector equation $Ax = \lambda x$.
 
-    By virtue of the Perron–Frobenius theorem [1]_, if G is strongly
-    connected there is a unique eigenvector $x$, and all its entries
+    By virtue of the Perron--Frobenius theorem [1]_, if `G` is (strongly)
+    connected, there is a unique eigenvector $x$, and all its entries
     are strictly positive.
 
-    If G is not strongly connected there might be several left
+    However, if `G` is not (strongly) connected, there might be several left
     eigenvectors associated with $\lambda$, and some of their elements
     might be zero.
+    Depending on the method used to choose eigenvectors, round-off error can affect
+    which of the infinitely many eigenvectors is reported.
+    This can lead to inconsistent results for the same graph,
+    which the underlying implementation is not robust to.
+    For this reason, only (strongly) connected graphs are accepted.
 
     Parameters
     ----------
     G : graph
-      A networkx graph.
+        A connected NetworkX graph.
+
+    weight : None or string, optional (default=None)
+        If ``None``, all edge weights are considered equal. Otherwise holds the
+        name of the edge attribute used as weight. In this measure the
+        weight is interpreted as the connection strength.
 
     max_iter : integer, optional (default=50)
-      Maximum number of Arnoldi update iterations allowed.
+        Maximum number of Arnoldi update iterations allowed.
 
     tol : float, optional (default=0)
-      Relative accuracy for eigenvalues (stopping criterion).
-      The default value of 0 implies machine precision.
-
-    weight : None or string, optional (default=None)
-      If None, all edge weights are considered equal. Otherwise holds the
-      name of the edge attribute used as weight. In this measure the
-      weight is interpreted as the connection strength.
+        Relative accuracy for eigenvalues (stopping criterion).
+        The default value of 0 implies machine precision.
 
     Returns
     -------
-    nodes : dictionary
-       Dictionary of nodes with eigenvector centrality as the value. The
-       associated vector has unit Euclidean norm and the values are
-       nonegative.
+    nodes : dict of nodes
+        Dictionary of nodes with eigenvector centrality as the value. The
+        associated vector has unit Euclidean norm and the values are
+        nonnegative.
 
     Examples
     --------
@@ -262,13 +267,16 @@ def eigenvector_centrality_numpy(G, weight=None, max_iter=50, tol=0):
     Raises
     ------
     NetworkXPointlessConcept
-        If the graph G is the null graph.
+        If the graph `G` is the null graph.
 
     ArpackNoConvergence
         When the requested convergence is not obtained. The currently
         converged eigenvalues and eigenvectors can be found as
         eigenvalues and eigenvectors attributes of the exception object.
 
+    AmbiguousSolution
+        If `G` is not connected.
+
     See Also
     --------
     :func:`scipy.sparse.linalg.eigs`
@@ -298,33 +306,35 @@ def eigenvector_centrality_numpy(G, weight=None, max_iter=50, tol=0):
     References
     ----------
     .. [1] Abraham Berman and Robert J. Plemmons.
-       "Nonnegative Matrices in the Mathematical Sciences."
+       "Nonnegative Matrices in the Mathematical Sciences".
        Classics in Applied Mathematics. SIAM, 1994.
 
     .. [2] Edmund Landau.
-       "Zur relativen Wertbemessung der Turnierresultate."
-       Deutsches Wochenschach, 11:366–369, 1895.
+       "Zur relativen Wertbemessung der Turnierresultate".
+       Deutsches Wochenschach, 11:366--369, 1895.
 
     .. [3] Teh-Hsing Wei.
-       "The Algebraic Foundations of Ranking Theory."
+       "The Algebraic Foundations of Ranking Theory".
        PhD thesis, University of Cambridge, 1952.
 
     .. [4] Maurice G. Kendall.
-       "Further contributions to the theory of paired comparisons."
-       Biometrics, 11(1):43–62, 1955.
+       "Further contributions to the theory of paired comparisons".
+       Biometrics, 11(1):43--62, 1955.
        https://www.jstor.org/stable/3001479
 
-    .. [5] Claude Berge
-       "Théorie des graphes et ses applications."
+    .. [5] Claude Berge.
+       "Théorie des graphes et ses applications".
        Dunod, Paris, France, 1958.
 
     .. [6] Phillip Bonacich.
-       "Technique for analyzing overlapping memberships."
-       Sociological Methodology, 4:176–185, 1972.
+       "Technique for analyzing overlapping memberships".
+       Sociological Methodology, 4:176--185, 1972.
        https://www.jstor.org/stable/270732
 
-    .. [7] Arnoldi iteration:: https://en.wikipedia.org/wiki/Arnoldi_iteration
-
+    .. [7] Arnoldi, W. E. (1951).
+       "The principle of minimized iterations in the solution of the matrix eigenvalue problem".
+       Quarterly of Applied Mathematics. 9 (1): 17--29.
+       https://doi.org/10.1090/qam/42792
     """
     import numpy as np
     import scipy as sp
@@ -333,6 +343,11 @@ def eigenvector_centrality_numpy(G, weight=None, max_iter=50, tol=0):
         raise nx.NetworkXPointlessConcept(
             "cannot compute centrality for the null graph"
         )
+    connected = nx.is_strongly_connected(G) if G.is_directed() else nx.is_connected(G)
+    if not connected:  # See gh-6888.
+        raise nx.AmbiguousSolution(
+            "`eigenvector_centrality_numpy` does not give consistent results for disconnected graphs"
+        )
     M = nx.to_scipy_sparse_array(G, nodelist=list(G), weight=weight, dtype=float)
     _, eigenvector = sp.sparse.linalg.eigs(
         M.T, k=1, which="LR", maxiter=max_iter, tol=tol
diff --git a/networkx/algorithms/centrality/reaching.py b/networkx/algorithms/centrality/reaching.py
index 93cb75a96..378e8a058 100644
--- a/networkx/algorithms/centrality/reaching.py
+++ b/networkx/algorithms/centrality/reaching.py
@@ -89,7 +89,6 @@ def global_reaching_centrality(G, weight=None, normalized=True):
     total_weight = G.size(weight=weight)
     if total_weight <= 0:
         raise nx.NetworkXError("Size of G must be positive")
-
     # If provided, weights must be interpreted as connection strength
     # (so higher weights are more likely to be chosen). However, the
     # shortest path algorithms in NetworkX assume the provided "weight"
@@ -178,10 +177,14 @@ def local_reaching_centrality(G, v, paths=None, weight=None, normalized=True):
            *PLoS ONE* 7.3 (2012): e33799.
            https://doi.org/10.1371/journal.pone.0033799
     """
+    # Corner case: graph with single node containing a self-loop
+    if (total_weight := G.size(weight=weight)) > 0 and len(G) == 1:
+        raise nx.NetworkXError(
+            "local_reaching_centrality of a single node with self-loop not well-defined"
+        )
     if paths is None:
         if nx.is_negatively_weighted(G, weight=weight):
             raise nx.NetworkXError("edge weights must be positive")
-        total_weight = G.size(weight=weight)
         if total_weight <= 0:
             raise nx.NetworkXError("Size of G must be positive")
         if weight is not None:
diff --git a/networkx/algorithms/distance_measures.py b/networkx/algorithms/distance_measures.py
index 0d8e52bf7..8e15bf8d9 100644
--- a/networkx/algorithms/distance_measures.py
+++ b/networkx/algorithms/distance_measures.py
@@ -1,11 +1,14 @@
 """Graph diameter, radius, eccentricity and other properties."""
 
+import math
+
 import networkx as nx
 from networkx.utils import not_implemented_for
 
 __all__ = [
     "eccentricity",
     "diameter",
+    "harmonic_diameter",
     "radius",
     "periphery",
     "center",
@@ -384,6 +387,72 @@ def diameter(G, e=None, usebounds=False, weight=None):
     return max(e.values())
 
 
+@nx._dispatchable
+def harmonic_diameter(G, sp=None):
+    """Returns the harmonic diameter of the graph G.
+
+    The harmonic diameter of a graph is the harmonic mean of the distances
+    between all pairs of distinct vertices. Graphs that are not strongly
+    connected have infinite diameter and mean distance, making such
+    measures not useful. Restricting the diameter or mean distance to
+    finite distances yields paradoxical values (e.g., a perfect match
+    would have diameter one). The harmonic mean handles gracefully
+    infinite distances (e.g., a perfect match has harmonic diameter equal
+    to the number of vertices minus one), making it possible to assign a
+    meaningful value to all graphs.
+
+    Note that in [1] the harmonic diameter is called "connectivity length":
+    however, "harmonic diameter" is a more standard name from the
+    theory of metric spaces. The name "harmonic mean distance" is perhaps
+    a more descriptive name, but is not used in the literature, so we use the
+    name "harmonic diameter" here.
+
+    Parameters
+    ----------
+    G : NetworkX graph
+       A graph
+
+    sp : dict of dicts, optional
+       All-pairs shortest path lengths as a dictionary of dictionaries
+
+    Returns
+    -------
+    hd : float
+       Harmonic diameter of graph
+
+    References
+    ----------
+    .. [1] Massimo Marchiori and Vito Latora, "Harmony in the small-world".
+           *Physica A: Statistical Mechanics and Its Applications*
+           285(3-4), pages 539-546, 2000.
+           <https://doi.org/10.1016/S0378-4371(00)00311-3>
+    """
+    order = G.order()
+
+    sum_invd = 0
+    for n in G:
+        if sp is None:
+            length = nx.single_source_shortest_path_length(G, n)
+        else:
+            try:
+                length = sp[n]
+                L = len(length)
+            except TypeError as err:
+                raise nx.NetworkXError('Format of "sp" is invalid.') from err
+
+        for d in length.values():
+            # Note that this will skip the zero distance from n to itself,
+            # as it should be, but also zero-weight paths in weighted graphs.
+            if d != 0:
+                sum_invd += 1 / d
+
+    if sum_invd != 0:
+        return order * (order - 1) / sum_invd
+    if order > 1:
+        return math.inf
+    return math.nan
+
+
 @nx._dispatchable(edge_attrs="weight")
 def periphery(G, e=None, usebounds=False, weight=None):
     """Returns the periphery of the graph G.


Test patch:
diff --git a/networkx/algorithms/centrality/tests/test_eigenvector_centrality.py b/networkx/algorithms/centrality/tests/test_eigenvector_centrality.py
index b8620056a..cfc9ee795 100644
--- a/networkx/algorithms/centrality/tests/test_eigenvector_centrality.py
+++ b/networkx/algorithms/centrality/tests/test_eigenvector_centrality.py
@@ -159,14 +159,26 @@ class TestEigenvectorCentralityExceptions:
         with pytest.raises(nx.NetworkXException):
             nx.eigenvector_centrality_numpy(nx.MultiGraph())
 
-    def test_empty(self):
+    def test_null(self):
         with pytest.raises(nx.NetworkXException):
             nx.eigenvector_centrality(nx.Graph())
 
-    def test_empty_numpy(self):
+    def test_null_numpy(self):
         with pytest.raises(nx.NetworkXException):
             nx.eigenvector_centrality_numpy(nx.Graph())
 
+    @pytest.mark.parametrize(
+        "G",
+        [
+            nx.empty_graph(3),
+            nx.DiGraph([(0, 1), (1, 2)]),
+        ],
+    )
+    def test_disconnected_numpy(self, G):
+        msg = "does not give consistent results for disconnected"
+        with pytest.raises(nx.AmbiguousSolution, match=msg):
+            nx.eigenvector_centrality_numpy(G)
+
     def test_zero_nstart(self):
         G = nx.Graph([(1, 2), (1, 3), (2, 3)])
         with pytest.raises(
diff --git a/networkx/algorithms/centrality/tests/test_reaching.py b/networkx/algorithms/centrality/tests/test_reaching.py
index 7267d8990..35d50e701 100644
--- a/networkx/algorithms/centrality/tests/test_reaching.py
+++ b/networkx/algorithms/centrality/tests/test_reaching.py
@@ -80,6 +80,17 @@ class TestGlobalReachingCentrality:
         actual = grc(G, normalized=False, weight="weight")
         assert expected == pytest.approx(actual, abs=1e-7)
 
+    def test_single_node_with_cycle(self):
+        G = nx.DiGraph([(1, 1)])
+        with pytest.raises(nx.NetworkXError, match="local_reaching_centrality"):
+            nx.global_reaching_centrality(G)
+
+    def test_single_node_with_weighted_cycle(self):
+        G = nx.DiGraph()
+        G.add_weighted_edges_from([(1, 1, 2)])
+        with pytest.raises(nx.NetworkXError, match="local_reaching_centrality"):
+            nx.global_reaching_centrality(G, weight="weight")
+
 
 class TestLocalReachingCentrality:
     """Unit tests for the local reaching centrality function."""
@@ -116,3 +127,14 @@ class TestLocalReachingCentrality:
             G, 1, normalized=True, weight="weight"
         )
         assert centrality == 1.0
+
+    def test_single_node_with_cycle(self):
+        G = nx.DiGraph([(1, 1)])
+        with pytest.raises(nx.NetworkXError, match="local_reaching_centrality"):
+            nx.local_reaching_centrality(G, 1)
+
+    def test_single_node_with_weighted_cycle(self):
+        G = nx.DiGraph()
+        G.add_weighted_edges_from([(1, 1, 2)])
+        with pytest.raises(nx.NetworkXError, match="local_reaching_centrality"):
+            nx.local_reaching_centrality(G, 1, weight="weight")
diff --git a/networkx/algorithms/tests/test_distance_measures.py b/networkx/algorithms/tests/test_distance_measures.py
index 97c3547f9..0b3840fd2 100644
--- a/networkx/algorithms/tests/test_distance_measures.py
+++ b/networkx/algorithms/tests/test_distance_measures.py
@@ -1,3 +1,4 @@
+import math
 from random import Random
 
 import pytest
@@ -52,6 +53,23 @@ class TestDistance:
     def test_diameter(self):
         assert nx.diameter(self.G) == 6
 
+    def test_harmonic_diameter(self):
+        assert abs(nx.harmonic_diameter(self.G) - 2.0477815699658715) < 1e-12
+
+    def test_harmonic_diameter_empty(self):
+        assert math.isnan(nx.harmonic_diameter(nx.empty_graph()))
+
+    def test_harmonic_diameter_single_node(self):
+        assert math.isnan(nx.harmonic_diameter(nx.empty_graph(1)))
+
+    def test_harmonic_diameter_discrete(self):
+        assert math.isinf(nx.harmonic_diameter(nx.empty_graph(3)))
+
+    def test_harmonic_diameter_not_strongly_connected(self):
+        DG = nx.DiGraph()
+        DG.add_edge(0, 1)
+        assert nx.harmonic_diameter(DG) == 2
+
     def test_radius(self):
         assert nx.radius(self.G) == 4