More efficient method for check_parallel_jacobian

idaes/core/util/model_diagnostics.py

@@ -27,7 +27,7 @@
 import numpy as np
 from scipy.linalg import svd
 from scipy.sparse.linalg import svds, norm
-from scipy.sparse import issparse, find
+from scipy.sparse import issparse, find, triu, diags
 
 from pyomo.environ import (
     Binary,
@@ -3619,7 +3619,109 @@ def ipopt_solve_halt_on_error(model, options=None):
     )
 
 
-def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "row"):
+def check_parallel_jacobian(
+    model,
+    tolerance: float = 1e-8,
+    direction: str = "row",
+    zero_norm_tolerance: float = 1e-8,
+    jac=None,
+    nlp=None,
+):
+    """
+    Check for near-parallel rows or columns in the Jacobian.
+
+    Near-parallel rows or columns indicate a potential degeneracy in the model,
+    as this means that the associated constraints or variables are (near)
+    duplicates of each other.
+
+    This method is based on work published in:
+
+    Klotz, E., Identification, Assessment, and Correction of Ill-Conditioning and
+    Numerical Instability in Linear and Integer Programs, Informs 2014, pgs. 54-108
+    https://pubsonline.informs.org/doi/epdf/10.1287/educ.2014.0130
+
+    Args:
+        model: model to be analysed
+        tolerance: tolerance to use to determine if constraints/variables are parallel
+        direction: 'row' (default, constraints) or 'column' (variables)
+
+    Returns:
+        list of 2-tuples containing parallel Pyomo components
+
+    """
+    # Thanks to Robby Parker for the sparse matrix implementation and
+    # significant performance improvements
+
+    if direction not in ["row", "column"]:
+        raise ValueError(
+            f"Unrecognised value for direction ({direction}). "
+            "Must be 'row' or 'column'."
+        )
+
+    if nlp is None or jac is None:
+        jac, nlp = get_jacobian(model, scaled=False)
+
+    # Get vectors that we will check, and the Pyomo components
+    # they correspond to.
+    if direction == "row":
+        components = nlp.get_pyomo_constraints()
+        mat = jac.tocsr()
+
+    elif direction == "column":
+        components = nlp.get_pyomo_variables()
+        mat = jac.transpose().tocsr()
+
+    # Take product of all rows/columns with all rows/columns by taking outer
+    # product of matrix with itself
+    outer = mat @ mat.transpose()
+
+    # Along the diagonal of the outer product you get the dot product of the row
+    # with itself, which is equal to the norm squared.
+    norms = np.sqrt(outer.diagonal())
+
+    zero_norm_indices = np.nonzero(np.abs(norms) <= zero_norm_tolerance)
+
+    # Want to ignore indices with zero norm. If we give them
+    # norms of infinity, we scale the rows and columns to zero
+    # which allows us to ignore them
+    norms[zero_norm_indices] = float("inf")  # 1/float('inf') == 0
+
+    # Divide each row and each column by the vector norm. This leaves
+    # the entries as dot(u, v) / (norm(u) * norm(v)). The exception is
+    # entries with "zero norm", whose corresponding rows and columns are
+    # set to zero.
+    scaling = diags(1 / norms)
+    outer = scaling * outer * scaling
+
+    # Get rid of duplicate values by only taking upper triangular part of
+    # resulting matrix
+    upper_tri = triu(outer)
+
+    # Set diagonal elements to zero
+    # Subtracting diags(upper_tri.diagonal()) is a more reliable
+    # method of getting the entries to exactly zero than subtracting
+    # an identity matrix, where one can encounter values of 1e-16
+    upper_tri = upper_tri - diags(upper_tri.diagonal())
+
+    # Get the nonzero entries of upper_tri in three arrays,
+    # corresponding to row indices, column indices, and values
+    rows, columns, values = find(upper_tri)
+
+    # We have that dot(u,v) == norm(u) * norm(v) * cos(theta) in which
+    # theta is the angle between u and v. If theta is approximately
+    # 0 or pi, sqrt(2*(1 - abs(dot(u,v)) / (norm(u) * norm(v)))) approximately
+    # equals the number of radians from 0 or pi. A tolerance of 1e-8 corresponds
+    # to a tolerance of sqrt(2) * 1e-4 radians
+    parallel_1D = np.nonzero(1 - abs(values) < tolerance)[0]
+
+    parallel = [
+        (components[rows[idx]], components[columns[idx]]) for idx in parallel_1D
+    ]
+
+    return parallel
+
+
+def check_parallel_jacobian_old(model, tolerance: float = 1e-4, direction: str = "row"):
     """
     Check for near-parallel rows or columns in the Jacobian.
 

idaes/core/util/tests/test_model_diagnostics.py

@@ -985,18 +985,12 @@ def test_display_near_parallel_variables(self):
 
         stream = StringIO()
         dt.display_near_parallel_variables(stream)
-
-        expected = """====================================================================================
-The following pairs of variables are nearly parallel:
-
-    v1, v2
-    v1, v4
-    v2, v4
-
-====================================================================================
-"""
-
-        assert stream.getvalue() == expected
+        expected_pairs = ["v1, v2", "v1, v4", "v1, v3", "v2, v4", "v2, v3", "v4, v3"]
+        assert (
+            "The following pairs of variables are nearly parallel:" in stream.getvalue()
+        )
+        for pair in expected_pairs:
+            assert pair in stream.getvalue()
 
     @pytest.mark.component
     def test_collect_structural_warnings_base_case(self, model):