diff --git a/idaes/core/util/linalg.py b/idaes/core/util/linalg.py
new file mode 100644
index 0000000000..b4358f4c57
--- /dev/null
+++ b/idaes/core/util/linalg.py
@@ -0,0 +1,392 @@
+#################################################################################
+# The Institute for the Design of Advanced Energy Systems Integrated Platform
+# Framework (IDAES IP) was produced under the DOE Institute for the
+# Design of Advanced Energy Systems (IDAES).
+#
+# Copyright (c) 2018-2025 by the software owners: The Regents of the
+# University of California, through Lawrence Berkeley National Laboratory,
+# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
+# University, West Virginia University Research Corporation, et al.
+# All rights reserved.  Please see the files COPYRIGHT.md and LICENSE.md
+# for full copyright and license information.
+#################################################################################
+"""
+This module contains specialized linear algebra routines not found in Numpy or Scipy
+"""
+
+__author__ = "Douglas Allan"
+
+import numpy as np
+from numpy.linalg import norm
+from numpy.random import default_rng
+from scipy.linalg import svd, eigh, qr
+from scipy.sparse.linalg import splu
+from scipy.sparse import issparse, block_array, eye as speye
+
+import idaes.logger as idaeslog
+
+_log = idaeslog.getLogger(__name__)
+
+
+def _symmetric_rayleigh_ritz_iteration(H, n_vec, tol, max_iter, seed=None):
+    """
+    Function to perform simultaneous Rayleigh quotient iteration on the
+    real symmetric matrix H.
+
+    Start out with initial eigenvalue estimates chosen randomly from
+    -1e-15 to 1e-15 in order to encourage convergence towards the
+    eigenvaleus of smallest magnitude, while the original eigenvector
+    estimates are chosen as an orthogonalization of a random matrix
+    whose values were chosen from the standard Gaussian distribution.
+
+    Then this iterative procedure is followed:
+        1) Compute the error in the eigenvalue approximations
+            err[j] = abs(H @ B[:, j] - mu[j] * B[:,j]) for all j
+        1a) If err[j] < tol for all j, terminate iteration
+        2) Choose j* such that err[j*] is maximized.
+        3) Calculate Bgrave = (H - mu[j*]) ** -1 @ B via a sparse
+            LU decomposition
+        4) Calculate an orthonormal basis Bhat for Bgrave using a
+            QR factorization
+        5) Calculate the Ritz values of H with respect to the
+            subspace Bhat via a symmetric eigendecomposition of
+            Bhat.T @ H @ Bhat and construct the Ritz vectors
+            from Bhat and the eigenvectors of Bhat.T @ H @ Bhat
+        6) Set mu and B to equal the Ritz values and vectors,
+            respectively, then return to (1).
+
+    The scalar version of this algorithm, in which B is a single
+    vector and mu is a scalar, is known as Rayleigh quotient
+    iteration, and is known to be cubically convergent for
+    symmetric matrices. I (Doug A.) have not seen this vectorized
+    extension anywhere in the literature. General attention
+    has been given to the (implicitly shifted) QR algorithm for
+    dense matrices and Krylov subspace methods for sparse matrices.
+    However, dense methods do not scale well and the Lanczos
+    algorithm (a Krylov method) provided by ARPACK through Scipy's
+    sparse.linalg.eigsh does not provide high quality eigenvector
+    approximations (or at least high enough quality approximations
+    to use as part of a sparse SVD routine), even when used in
+    shift-invert mode on the augmented matrix.
+
+    For finding small eigenvalues, the methods provided in Section 5
+    of Sleijpen and Van Der Vorst (2000) and Section 4.4 of
+    Hochstenbach (2001) may provide performant alternatives, because
+    they can be implemented in a vectorized form in Python. The
+    Lanczos variant suggested by Kokiopoulou et al. (2004) may
+    be worth attention, but may have to be implemented in a lower
+    level language to be performant.
+
+    Sleijpen, G.L.G., Van Der Vorst, H.A., 2000. A Jacobi--Davidson
+    Iteration Method for Linear Eigenvalue Problems. SIAM Rev. 42,
+    267–293. https://doi.org/10.1137/S0036144599363084
+
+    Hochstenbach, M.E., 2001. A Jacobi--Davidson Type SVD Method.
+    SIAM J. Sci. Comput. 23, 606–628.
+    https://doi.org/10.1137/S1064827500372973
+
+    Kokiopoulou, E., Bekas, C., Gallopoulos, E., 2004. Computing
+    smallest singular triplets with implicitly restarted Lanczos
+    bidiagonalization. Applied Numerical Mathematics 49, 39–61.
+    https://doi.org/10.1016/j.apnum.2003.11.011
+
+
+
+    Args:
+        H: Real, symmetric n by n matrix H
+        n_vec: Number of eigenvectors and eigenvalues to calculate
+        tol: Tolerance used to decide to terminate iteration.
+        max_iter: maximum number of iterations
+
+    Returns:
+        evals: 1D array of n_vec eigenvalues
+        evecs: n by n_vec 2D array of eigenvectors
+        converged: Boolean flag of whether the convergence criteria
+            was met after max_iter iterations
+    """
+    assert len(H.shape) == 2
+    m, n = H.shape
+    assert m == n
+    rng_obj = default_rng(seed)
+    mu = rng_obj.uniform(low=-1e-15, high=1e-15, size=(n_vec,))
+    B = rng_obj.standard_normal(size=(m, n_vec))
+    B, _ = qr(B, mode="economic")
+    # import pdb; pdb.set_trace()
+
+    converged = False
+
+    for i in range(max_iter):
+        # Calculate how close B[:,k] and mu[k] are to an eigenpair
+        # by calculating norm(H @ B[:, k]  - mu[k] * B[:,k]) for all k.
+        err = norm(H @ B - B @ np.diag(mu), axis=0)
+
+        # TODO should this tolerance be scaled by matrix size?
+        if max(err) < tol:
+            converged = True
+            print(f"Converged in {i} iterations")
+            break
+
+        # Find the index which is furthest from being an eigenvector in
+        # order to improve the estimate for that index
+        max_err_idx = np.argmax(err)
+
+        try:
+            # Calculate
+            invH_shift = splu(H - mu[max_err_idx] * speye(n, format="csc"))
+        except RuntimeError as exc1:
+            if "Factor is exactly singular" in str(exc1):
+                # mu[max_err_idx] is exactly equal to an
+                # eigenvalue, but B[:,max_err_idx] is not
+                # a good eigenvector estimate. Therefore
+                # perturb mu[max_err_idx] to allow matrix
+                # inversion to succeed
+                try:
+                    invH_shift = splu(
+                        H - (1e-15 + mu[max_err_idx]) * speye(n, format="csc")
+                    )
+                except RuntimeError as exc2:
+                    if "Factor is exactly singular" in str(exc2):
+                        # If perturbation fails, we'll assume
+                        # something bigger is wrong in the solver
+                        raise RuntimeError(
+                            "Unable to invert matrix H when shifted by "
+                            f"{mu[max_err_idx]}*I. Check whether something "
+                            "is wrong with the matrix H's scaling."
+                        ) from exc2
+                    else:
+                        raise
+            else:
+                raise
+        Bgrave = invH_shift.solve(B)
+        Bhat, _ = qr(Bgrave, mode="economic")
+        mu, B_tilde = eigh(Bhat.T @ H @ Bhat)
+        B = Bhat @ B_tilde
+
+    return mu, B, converged
+
+
+def _aug_eig_processing(
+    A,
+    evecs,
+    zero_tol=1e-2,
+):
+    """
+    Takes the output of subspace iteration on the augmented matrix [[0, A.T], [A, 0]]
+    for some m by n matrix A, and obtains estimates for U, svals, V such that
+    U.T @ A @ V == diag(svals). If m < n (m > n), then a basis for the (left) null
+    space is also computed.
+
+    For every singular triplet (sigma, u, v), the augmented matrix has two eigentriplets:
+    (sigma, (v, u)) and (-sigma, (-v, u)). If m!=n, there are |m-n| eigenvectors corresponding
+    to the null space. If m > n, then these vectors are of the form (0, (0, w)), and if m < n
+    they are of the form (0, (w, 0)). We take eigenvectors of the augmented matrix as an input,
+    but there is no telling how many singular vectors they correspond to: sometimes the method
+    converges to both (sigma, (v, u)) and (-sigma, (-v, u)) and sometimes it only converges to
+    a single one. Furthermore, if some sigma is close to 0, we can end up with blending between
+    the null space and the space corresponding to the small singular value. Finally, while the
+    eigenvectors may be orthogonal to machine precision, the resulting singular vectors may not
+    be orthogonal due to error cancellation between the u and v components: for (v_1, u_1) and
+    (v_2, u_2), we may have v_1.T @ v_2 + u_1.T @ u_2 ~= 1e-14, but nevertheless have
+    v_1.T @ v_2 ~= 1e-7 and u_1.T @ u_2 ~= 1e-7.
+
+    The solution to this very messy situation is, as usual, the QR factorization. We partition
+    the eigenvectors into U and V components, then perform a rank-revealing QR in order to
+    obtain orthonormal bases of those sets of vectors. The diagonal elements of R are measures
+    of how strongly a basis vector is represented in the U and V components of evecs. In an
+    ideal world, if evecs contained q singular vectors and p null vectors, we would have
+    R[k,k] == 1 for k<=p-1 (corresponding to the null vectors), R[k,k] = 1/sqrt(2) for
+    p <= k <= p+q-1 (corresponding to the singular vectors), and R[k,k] ~= 1e-15 for k > p+q.
+    However, we often get a spectrum of values.
+
+    Presently, we use a crude thresholding scheme: if both abs(R_U[k,k]) and abs(R_V[k,k]) are
+    greater than zero_tol, we add the columns Q_U[:,k] and Q_V[:,k] to our basis. If only one
+    is greater than zero_tol, then we check how many more basis vectors we have for U than for
+    V (if m > n) or for V than U (if m < n). If the difference in number of basis vectors is
+    less than the size of the null space, we add an additional basis vector, otherwise we drop
+    it. This scheme is somewhat arbitrary, but has worked better than any replacement that I
+    (Doug A.) have tried.
+
+    After determining bases for U and V, we perform a SVD on U.T @ A @ V in order to
+    reconstruct appropriate singular triplets as well as the null vectors.
+
+    Args:
+        A: Scipy sparse array with real entries
+        evecs: Approximate eigenvectors of the augmented matrix [[0, A.T], [A, 0]]
+        zero_tol: Tolerance to drop vectors for basis of singular spaces U and V
+
+    Returns:
+        U: Dense array of left singular vectors of A
+        svals: 1D dense array ofsingular values of A
+        V: Dense array of left singular vectors of A
+        null_space: If m < n (m > n) then a basis for the (left) null space
+            is returned. If m == n, then no null space is returned.
+
+    """
+    assert len(A.shape) == 2
+    m, n = A.shape
+    p = abs(m - n)
+
+    V = evecs[:n, :]
+    U = evecs[n:, :]
+
+    U, R_U, _ = qr(U, mode="economic", pivoting=True)
+    r_U = np.abs(np.diag(R_U))
+
+    V, R_V, _ = qr(V, mode="economic", pivoting=True)
+    r_V = np.abs(np.diag(R_V))
+
+    n_U = 0
+    n_V = 0
+    for k in range(min(evecs.shape[1], max(m, n))):
+        if m > n:
+            if r_U[k] > zero_tol:
+                if k < n and r_V[k] > zero_tol:
+                    n_U += 1
+                    n_V += 1
+                elif n_U - n_V < p:
+                    n_U += 1
+        else:
+            if r_V[k] > zero_tol:
+                if k < m and r_U[k] > zero_tol:
+                    n_U += 1
+                    n_V += 1
+
+                elif n_V - n_U < p:
+                    n_V += 1
+
+    U = U[:, :n_U]
+    V = V[:, :n_V]
+
+    # Want to make sure that the number of null vector candidates
+    # is the expected size
+    # This assertion should no longer be necessary
+    p_hat = abs(U.shape[1] - V.shape[1])
+    assert p_hat <= p
+    U_sub, svals, VT_sub = svd(U.T @ A @ V, full_matrices=True)
+    U = U @ U_sub
+    V = V @ VT_sub.T
+
+    if m > n:
+        null_space = U[:, -p_hat:]
+        assert norm(null_space.T @ A) <= np.sqrt(m * n) * 1e-15
+        U = U[:, :-p_hat]
+    elif m < n:
+        null_space = V[:, -p_hat:]
+        assert norm(A @ null_space) <= np.sqrt(m * n) * 1e-15
+        V = V[:, :-p_hat]
+
+    # Sort singular values in ascending order
+    U = U[:, ::-1]
+    V = V[:, ::-1]
+    svals = svals[::-1]
+
+    if m == n:
+        return U, svals, V
+    else:
+        return U, svals, V, null_space
+
+
+def svd_rayleigh_ritz(
+    A,
+    number_singular_values: int = 10,
+    max_iter: int = 100,
+    tol: float = 1e-12,
+    seed: int = None,
+    suppress_warning=False,
+):
+    """
+    Computes smallest singular vectors of the sparse m by n matrix A via Rayleigh-Ritz
+    iteration on the augmented matrix [[0, A.T], [A, 0]]. Working with this matrix
+    avoids the roundoff error inherent working with the Gram matrices A.T @ A or
+    A @ A.T, but also results in the (left) null space polluting the singular spectrum
+    if m < n (m > n). Therefore, this method is not appropriate for diagnosing
+    optimization problems in which m << n, but is appropriate if the (left) null space
+    is desired for diagnosing degrees of freedom.
+
+    Args:
+        A: Scipy sparse array with real entries
+        number_singular_values: number of small singular values and vectors to compute
+        max_iter: Maximum number of iterations for Rayleigh-Ritz iteration
+        tol: Tolerance used in stopping condition for Rayleigh-Ritz iteration
+        seed: Seed for initializing random number generator in Rayleigh-Ritz iteration
+        suppress_warning: Suppress the efficiency warning issued when |m - n| > 10
+
+    Returns:
+        U: m by number_singular_values dense array of left singular vectors
+        svals: 1D dense array of number_singular_values singular values
+        V: n by number_singular_values dense array of left singular vectors
+        null: Basis for the (left) null space if m < n (m > n)
+    """
+    # Should we try to squeeze extra dimensions first? It doesn't look like np.squeeze
+    # works on sparse arrays, so we'd need to reshape it.
+
+    # It appears that my version of Scipy does not support sparse arrays of dimension
+    # higher than 2, but that such a feature is actively being worked on. Therefore
+    # the shape of A first to facilitate testing
+    if not len(A.shape) == 2:
+        raise ValueError(
+            "This method expects a 2D Scipy sparse array-like as input, but was passed "
+            f"a {len(A.shape)}D array-like instead."
+        )
+
+    if not issparse(A):
+        raise ValueError(
+            "This method expects a Scipy sparse array-like as an input but was passed "
+            "a dense array-like instead. Try using scipy.linalg.svd for a dense SVD method."
+        )
+    m, n = A.shape
+    l = abs(m - n)
+    if l > 10 and not suppress_warning:
+        if m > n:
+            _log.warning(
+                f"Matrix A has a left nullspace of dimension at least {l}, which "
+                "degrades the efficiency of SVD algorithms based on the augmented "
+                "matrix [[0, A.T], [A, 0]]."
+            )
+        else:
+            _log.warning(
+                f"Matrix A has a nullspace of dimension at least {l}, which "
+                "degrades the efficiency of SVD algorithms based on the augmented "
+                "matrix [[0, A.T], [A, 0]]."
+            )
+    # Can't obtain more singular values than the dimension of the matrix
+    number_singular_values = min(number_singular_values, m, n)
+
+    # The augmented matrix has two eigenvectors for every singular triplet, in addition
+    # to eigenvectors corresponding to the (left) null space if m<n (m>n). In order to
+    # get the required number of singular triplets, we need to exhaust the null space
+    # as well as the duplicate singular triplets.
+    n_samples = 2 * number_singular_values + l
+    A_aug = block_array([[None, A.T], [A, None]], format="csc")
+
+    _, B, converged = _symmetric_rayleigh_ritz_iteration(
+        A_aug, n_samples, tol=tol, max_iter=max_iter, seed=seed
+    )
+
+    if not converged:
+        raise RuntimeError(
+            "Rayleigh-Ritz iteration did not converge! Consider increasing the tolerance or "
+            "maximum number of iterations."
+        )
+
+    try:
+        if m == n:
+            U, svals, V = _aug_eig_processing(A, B)
+        else:
+            U, svals, V, null = _aug_eig_processing(A, B)
+    except AssertionError as exc:
+        raise RuntimeError(
+            "Processing of singular vectors failed despite Rayleigh-Ritz iteration "
+            "converging. Please let the IDAES dev team know about this failure so that "
+            "this processing step can be refined."
+        ) from exc
+
+    # Singular values already in ascending order, so just take the number that we want
+    U = U[:, :number_singular_values]
+    V = V[:, :number_singular_values]
+    svals = svals[:number_singular_values]
+
+    if m == n:
+        return U, svals, V
+    else:
+        return U, svals, V, null
diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py
index fe5af973ee..d2821e102e 100644
--- a/idaes/core/util/model_diagnostics.py
+++ b/idaes/core/util/model_diagnostics.py
@@ -124,6 +124,8 @@
     ParameterSweepBase,
     is_psweepspec,
 )
+from idaes.core.util.linalg import svd_rayleigh_ritz
+
 import idaes.logger as idaeslog
 
 _log = idaeslog.getLogger(__name__)
@@ -200,6 +202,29 @@ def svd_sparse(jacobian, number_singular_values):
     return u, s, vT.transpose()
 
 
+def svd_rayleigh_ritz_callback(jacobian, number_singular_values, **kwargs):
+    """
+    Callback for performing SVD analysis using idaes.core.util.linalg.svd_rayleigh_ritz
+
+    Args:
+        jacobian: Jacobian to be analysed
+        number_singular_values: number of singular values to compute
+        **kwargs: Dictionary of keyword arguments to pass to svd_rayleigh_ritz
+
+    Returns:
+        u, s and v numpy arrays
+
+    """
+    # This method also returns the null space, which is not used by
+    # the model diagnostics at present
+    m, n = jacobian.shape
+    if m != n:
+        u, s, v, _ = svd_rayleigh_ritz(jacobian, number_singular_values, **kwargs)
+    else:
+        u, s, v = svd_rayleigh_ritz(jacobian, number_singular_values, **kwargs)
+    return u, s, v
+
+
 CONFIG = ConfigDict()
 CONFIG.declare(
     "variable_bounds_absolute_tolerance",
@@ -363,10 +388,10 @@ def svd_sparse(jacobian, number_singular_values):
 SVDCONFIG.declare(
     "svd_callback",
     ConfigValue(
-        default=svd_dense,
+        default=svd_rayleigh_ritz_callback,
         domain=svd_callback_validator,
-        description="Callback to SVD method of choice (default = svd_dense)",
-        doc="Callback to SVD method of choice (default = svd_dense). "
+        description="Callback to SVD method of choice (default = svd_rayleigh_ritz_callback)",
+        doc="Callback to SVD method of choice (default = svd_rayleigh_ritz_callback). "
         "Callbacks should take the Jacobian and number of singular values "
         "to compute as options, plus any method specific arguments, and should "
         "return the u, s and v matrices as numpy arrays.",
diff --git a/idaes/core/util/tests/svd_cache/overdetermined_mea_column_jac.npz b/idaes/core/util/tests/svd_cache/overdetermined_mea_column_jac.npz
new file mode 100644
index 0000000000..398dd6be36
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/overdetermined_mea_column_jac.npz differ
diff --git a/idaes/core/util/tests/svd_cache/overdetermined_mea_column_svd.npz b/idaes/core/util/tests/svd_cache/overdetermined_mea_column_svd.npz
new file mode 100644
index 0000000000..18cff96fc8
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/overdetermined_mea_column_svd.npz differ
diff --git a/idaes/core/util/tests/svd_cache/overdetermined_soc_jac.npz b/idaes/core/util/tests/svd_cache/overdetermined_soc_jac.npz
new file mode 100644
index 0000000000..08d82cca99
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/overdetermined_soc_jac.npz differ
diff --git a/idaes/core/util/tests/svd_cache/overdetermined_soc_svd.npz b/idaes/core/util/tests/svd_cache/overdetermined_soc_svd.npz
new file mode 100644
index 0000000000..26c7c4633a
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/overdetermined_soc_svd.npz differ
diff --git a/idaes/core/util/tests/svd_cache/square_big_mea_column_jac.npz b/idaes/core/util/tests/svd_cache/square_big_mea_column_jac.npz
new file mode 100644
index 0000000000..c92a5f326a
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/square_big_mea_column_jac.npz differ
diff --git a/idaes/core/util/tests/svd_cache/square_big_mea_column_svd.npz b/idaes/core/util/tests/svd_cache/square_big_mea_column_svd.npz
new file mode 100644
index 0000000000..4fe115bca2
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/square_big_mea_column_svd.npz differ
diff --git a/idaes/core/util/tests/svd_cache/underdetermined_mea_column_jac.npz b/idaes/core/util/tests/svd_cache/underdetermined_mea_column_jac.npz
new file mode 100644
index 0000000000..1cd54c2963
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/underdetermined_mea_column_jac.npz differ
diff --git a/idaes/core/util/tests/svd_cache/underdetermined_mea_column_svd.npz b/idaes/core/util/tests/svd_cache/underdetermined_mea_column_svd.npz
new file mode 100644
index 0000000000..83365e5a31
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/underdetermined_mea_column_svd.npz differ
diff --git a/idaes/core/util/tests/svd_cache/underdetermined_soc_jac.npz b/idaes/core/util/tests/svd_cache/underdetermined_soc_jac.npz
new file mode 100644
index 0000000000..a856941337
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/underdetermined_soc_jac.npz differ
diff --git a/idaes/core/util/tests/svd_cache/underdetermined_soc_svd.npz b/idaes/core/util/tests/svd_cache/underdetermined_soc_svd.npz
new file mode 100644
index 0000000000..79e2790985
Binary files /dev/null and b/idaes/core/util/tests/svd_cache/underdetermined_soc_svd.npz differ
diff --git a/idaes/core/util/tests/test_linalg.py b/idaes/core/util/tests/test_linalg.py
new file mode 100644
index 0000000000..d007bdde34
--- /dev/null
+++ b/idaes/core/util/tests/test_linalg.py
@@ -0,0 +1,372 @@
+#################################################################################
+# The Institute for the Design of Advanced Energy Systems Integrated Platform
+# Framework (IDAES IP) was produced under the DOE Institute for the
+# Design of Advanced Energy Systems (IDAES).
+#
+# Copyright (c) 2018-2024 by the software owners: The Regents of the
+# University of California, through Lawrence Berkeley National Laboratory,
+# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
+# University, West Virginia University Research Corporation, et al.
+# All rights reserved.  Please see the files COPYRIGHT.md and LICENSE.md
+# for full copyright and license information.
+#################################################################################
+"""
+Tests for linear algebra utility methods.
+"""
+
+import os
+import pytest
+
+import numpy as np
+from numpy.linalg import norm, qr
+from numpy.random import default_rng
+from scipy.sparse import csc_array, load_npz
+
+from pyomo.common.fileutils import this_file_dir
+from pyomo.environ import log10
+
+import idaes.logger as idaeslog
+from idaes.core.util.linalg import svd_rayleigh_ritz
+
+svd_cache = os.sep.join([this_file_dir(), "svd_cache"])
+
+__author__ = "Douglas Allan"
+
+
+def _random_svd(
+    n_rows: int,
+    n_columns: int,
+    n_small: int,
+    seed: int,
+    eps_small: float = 1e-12,
+    eps_min: float = 1e-16,
+):
+    """
+    Constructs an n_rows by n_columns Numpy 2D array A from its singular value decomposition
+    A = U @ diag(sigma) @ V.T. U and V are created through orthogonalization of random
+    arrays. The 1D array sigma is generated such that min(n_rows, n_columns) - n_small
+    singular values are greater than eps_small while n_small singular values are less than
+    eps_small.
+
+    Args:
+        n_rows: Number of rows in outlet array A
+        n_columns: Number of columns in outlet array A
+        n_small: Number of singular values smaller than eps_small
+        seed: Seed for Numpy's random number generator
+        eps_small: Value defining threshold for small singular value (default 1e-12)
+        eps_min: Smallest singular value possible
+    Returns:
+        A: Random matrix generated by function
+        U: m by number_singular_values dense array of left singular vectors
+        svals: 1D dense array of number_singular_values singular values, sorted from smallest to largest
+        V: n by number_singular_values dense array of left singular vectors
+    """
+    n_svals = min(n_rows, n_columns)
+    assert n_svals >= n_small
+    rng_obj = default_rng(seed)
+
+    U = rng_obj.standard_normal((n_rows, n_svals))
+    V = rng_obj.standard_normal((n_columns, n_svals))
+    U, _ = qr(U)
+    V, _ = qr(V)
+
+    svals_big = rng_obj.uniform(low=log10(eps_small), high=0, size=(n_svals - n_small,))
+    svals_big = 10**svals_big
+    svals_big.sort()
+    svals_small = rng_obj.uniform(
+        low=log10(eps_min), high=log10(eps_small), size=(n_small,)
+    )
+    svals_small = 10**svals_small
+    svals_small.sort()
+    svals = np.concatenate([svals_small, svals_big])
+
+    A = U @ np.diag(svals) @ V.T
+
+    return A, U, svals, V
+
+
+def _assert_subspace_containment(U, V, tol=1e-8):
+    """
+    Asserts that the matrix U is contained in the subspace spanned by
+    the orthonormal basis V to an absolute tolerance of tol.
+    """
+    assert norm(V @ (V.T @ U) - U) == pytest.approx(0, rel=0, abs=tol)
+
+
+def _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat=None, nvec=10):
+    m, n = A.shape
+    abs_tol = 1e-15 * np.sqrt(m * n)
+    # First test shapes
+    assert Uhat.shape == (m, nvec)
+    assert Vhat.shape == (n, nvec)
+    assert svals_hat.shape == (nvec,)
+    if null_hat is not None:
+        # It's highly unlikely, but we may get fewer null vectors than |m - n|
+        # so we don't test for that here.
+        if m > n:
+            assert null_hat.shape[0] == m
+        else:
+            assert null_hat.shape[0] == n
+
+    # U.T @ A @ V = Sigma
+    assert Uhat.T @ A @ Vhat == pytest.approx(np.diag(svals_hat), rel=1e-6, abs=abs_tol)
+    # U.T @ A = Sigma @ V.T
+    assert norm(Uhat.T @ A, axis=1) == pytest.approx(svals_hat, rel=1e-2, abs=abs_tol)
+    # A @ V = U @ Sigma
+    assert norm(A @ Vhat, axis=0) == pytest.approx(svals_hat, rel=1e-2, abs=abs_tol)
+    if null_hat is not None:
+        if m > n:
+            assert norm(null_hat.T @ A) == pytest.approx(0, rel=0, abs=abs_tol)
+        else:
+            assert norm(A @ null_hat) == pytest.approx(0, rel=0, abs=abs_tol)
+
+    # Check orthogonality
+    assert Uhat.T @ Uhat == pytest.approx(np.eye(Uhat.shape[1]), rel=0, abs=abs_tol)
+    assert Vhat.T @ Vhat == pytest.approx(np.eye(Vhat.shape[1]), rel=0, abs=abs_tol)
+    if null_hat is not None:
+        assert null_hat.T @ null_hat == pytest.approx(
+            np.eye(null_hat.shape[1]), rel=0, abs=abs_tol
+        )
+        if m > n:
+            assert norm(null_hat.T @ Uhat) == pytest.approx(0, abs=abs_tol)
+        else:
+            assert norm(null_hat.T @ Vhat) == pytest.approx(0, abs=abs_tol)
+
+
+class TestSVDRayleighRitz:
+    @pytest.mark.unit
+    def test_not_sparse(self):
+        A = np.eye(2)
+        with pytest.raises(
+            ValueError,
+            match="This method expects a Scipy sparse array-like as an input but was passed "
+            "a dense array-like instead. Try using scipy.linalg.svd for a dense SVD method.",
+        ):
+            svd_rayleigh_ritz(A)
+
+    @pytest.mark.unit
+    def test_not_2D(self):
+        A = np.zeros((2, 3, 4))
+        with pytest.raises(
+            ValueError,
+            match="This method expects a 2D Scipy sparse array-like as input, but was passed "
+            "a 3D array-like instead.",
+        ):
+            svd_rayleigh_ritz(A)
+
+    @pytest.mark.unit
+    def test_not_converged(self):
+        A = np.eye(2)
+        with pytest.raises(
+            RuntimeError,
+            match="Rayleigh-Ritz iteration did not converge! Consider increasing "
+            "the tolerance or maximum number of iterations.",
+        ):
+            svd_rayleigh_ritz(csc_array(A), max_iter=0)
+
+    @pytest.mark.unit
+    def test_square(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=30, n_columns=30, n_small=5, seed=7, eps_min=1e-16
+        )
+        svals_trunc = svals[:10]
+        Uhat, svals_hat, Vhat = svd_rayleigh_ritz(csc_array(A), seed=1)
+        _test_svd_quality(A, Uhat, svals_hat, Vhat)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_overdetermined(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=35, n_columns=30, n_small=5, seed=8, eps_min=1e-16
+        )
+        svals_trunc = svals[:10]
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(csc_array(A), seed=2)
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_underdetermined(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=30, n_columns=35, n_small=5, seed=9, eps_min=1e-16
+        )
+        svals_trunc = svals[:10]
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(csc_array(A), seed=3)
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_underdetermined_warning(self, caplog):
+        A, _, svals, _ = _random_svd(
+            n_rows=30, n_columns=45, n_small=5, seed=13, eps_min=1e-16
+        )
+        svals_trunc = svals[:10]
+        with caplog.at_level(idaeslog.WARNING):
+            Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(csc_array(A), seed=3)
+        assert (
+            "Matrix A has a nullspace of dimension at least 15, which "
+            "degrades the efficiency of SVD algorithms based on the augmented "
+            "matrix [[0, A.T], [A, 0]]."
+        ) in caplog.text
+
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_overdetermined_warning(self, caplog):
+        A, _, svals, _ = _random_svd(
+            n_rows=44, n_columns=30, n_small=5, seed=23, eps_min=1e-16
+        )
+        svals_trunc = svals[:10]
+        with caplog.at_level(idaeslog.WARNING):
+            Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(csc_array(A), seed=3)
+        assert (
+            "Matrix A has a left nullspace of dimension at least 14, which "
+            "degrades the efficiency of SVD algorithms based on the augmented "
+            "matrix [[0, A.T], [A, 0]]."
+        ) in caplog.text
+
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_underdetermined_small(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=5, n_columns=6, n_small=1, seed=657457, eps_min=1e-16
+        )
+        svals_trunc = svals
+        # The default option is 10 singular vectors, but the method is supposed to silently
+        # truncate it to 5 singular vectors
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(csc_array(A), seed=12342314)
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat, nvec=5)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_overdetermined_small(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=6, n_columns=5, n_small=1, seed=4321412, eps_min=1e-16
+        )
+        svals_trunc = svals
+        # The default option is 10 singular vectors, but the method is supposed to silently
+        # truncate it to 5 singular vectors
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(csc_array(A), seed=3425767)
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat, nvec=5)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_square_small(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=5, n_columns=5, n_small=1, seed=5463, eps_min=1e-16
+        )
+        svals_trunc = svals
+        # The default option is 10 singular vectors, but the method is supposed to silently
+        # truncate it to 5 singular vectors
+        Uhat, svals_hat, Vhat = svd_rayleigh_ritz(csc_array(A), seed=9876575)
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, nvec=5)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_more_svs(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=30, n_columns=35, n_small=5, seed=12, eps_min=1e-16
+        )
+        svals_trunc = svals[:12]
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(
+            csc_array(A), number_singular_values=12, seed=31
+        )
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat, nvec=12)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.unit
+    def test_fewer_svs(self):
+        A, _, svals, _ = _random_svd(
+            n_rows=30, n_columns=35, n_small=5, seed=37, eps_min=1e-16
+        )
+        svals_trunc = svals[:7]
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(
+            csc_array(A), number_singular_values=7, seed=31
+        )
+        _test_svd_quality(A, Uhat, svals_hat, Vhat, null_hat, nvec=7)
+        assert pytest.approx(svals_trunc, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.component
+    def test_soc_overdetermined(self):
+        jac = load_npz(os.sep.join([svd_cache, "overdetermined_soc_jac.npz"]))
+        # "Truth" values come from a dense svd of the matrix
+        cached_svd = np.load(
+            os.sep.join([svd_cache, "overdetermined_soc_svd.npz"]), allow_pickle=False
+        )
+        Utrue = cached_svd["Utrue"]
+        svals_true = cached_svd["svals_true"]
+        Vtrue = cached_svd["Vtrue"]
+        null_true = cached_svd["null_true"]
+
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(jac, seed=86)
+        _test_svd_quality(jac, Uhat, svals_hat, Vhat, null_hat)
+        _assert_subspace_containment(Uhat, Utrue)
+        _assert_subspace_containment(Vhat, Vtrue)
+        _assert_subspace_containment(null_hat, null_true)
+        assert pytest.approx(svals_true, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.integration
+    def test_soc_underdetermined(self):
+        jac = load_npz(os.sep.join([svd_cache, "underdetermined_soc_jac.npz"]))
+        # "Truth" values come from a dense svd of the matrix
+        cached_svd = np.load(
+            os.sep.join([svd_cache, "underdetermined_soc_svd.npz"]), allow_pickle=False
+        )
+        Utrue = cached_svd["Utrue"]
+        svals_true = cached_svd["svals_true"]
+        Vtrue = cached_svd["Vtrue"]
+        null_true = cached_svd["null_true"]
+
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(jac, seed=86, max_iter=200)
+        _test_svd_quality(jac, Uhat, svals_hat, Vhat, null_hat)
+        _assert_subspace_containment(Uhat, Utrue)
+        _assert_subspace_containment(Vhat, Vtrue)
+        _assert_subspace_containment(null_hat, null_true)
+        assert pytest.approx(svals_true, rel=1e-4, abs=1e-14) == svals_hat
+
+    @pytest.mark.integration
+    def test_mea_column_overdetermined(self):
+        jac = load_npz(os.sep.join([svd_cache, "overdetermined_mea_column_jac.npz"]))
+        # "Truth" values come from a dense svd of the matrix
+        # However, with this seed, we do not converge to the 10 smallest singular
+        # values, but only 9 of them plus the 11th(?) smallest sval
+
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(jac, seed=42, max_iter=200)
+        _test_svd_quality(jac, Uhat, svals_hat, Vhat, null_hat)
+
+    @pytest.mark.integration
+    def test_mea_column_underdetermined(self):
+        jac = load_npz(os.sep.join([svd_cache, "underdetermined_mea_column_jac.npz"]))
+        # "Truth" values come from a dense svd of the matrix
+        # We've got the smallest 20 singular vectors saved, so we can test subspace containment here
+        cached_svd = np.load(
+            os.sep.join([svd_cache, "underdetermined_mea_column_svd.npz"]),
+            allow_pickle=False,
+        )
+        Utrue = cached_svd["Utrue"]
+        Vtrue = cached_svd["Vtrue"]
+        null_true = cached_svd["null_true"]
+
+        Uhat, svals_hat, Vhat, null_hat = svd_rayleigh_ritz(jac, seed=87, max_iter=200)
+        _test_svd_quality(jac, Uhat, svals_hat, Vhat, null_hat)
+        _assert_subspace_containment(Uhat, Utrue)
+        _assert_subspace_containment(Vhat, Vtrue)
+        _assert_subspace_containment(null_hat, null_true)
+
+    @pytest.mark.integration
+    def test_big_mea_column_square(self):
+        jac = load_npz(os.sep.join([svd_cache, "square_big_mea_column_jac.npz"]))
+        # "Truth" values come from a dense svd of the matrix
+        cached_svd = np.load(
+            os.sep.join([svd_cache, "square_big_mea_column_svd.npz"]),
+            allow_pickle=False,
+        )
+        Utrue = cached_svd["Utrue"]
+        Vtrue = cached_svd["Vtrue"]
+
+        Uhat, svals_hat, Vhat = svd_rayleigh_ritz(jac, seed=87, max_iter=200)
+        _test_svd_quality(jac, Uhat, svals_hat, Vhat)
+        _assert_subspace_containment(Uhat, Utrue, tol=1e-6)
+        _assert_subspace_containment(Vhat, Vtrue, tol=1e-6)
diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py
index 4f402dc73b..5029b1f1a0 100644
--- a/idaes/core/util/tests/test_model_diagnostics.py
+++ b/idaes/core/util/tests/test_model_diagnostics.py
@@ -87,6 +87,7 @@
     DegeneracyHunter2,
     svd_dense,
     svd_sparse,
+    svd_rayleigh_ritz_callback,
     get_valid_range_of_component,
     set_bounds_from_valid_range,
     list_components_with_values_outside_valid_range,
@@ -2119,9 +2120,29 @@ def test_init_small_model(self):
     @pytest.mark.unit
     def test_run_svd_analysis(self, dummy_problem):
         svd = SVDToolbox(dummy_problem)
+        assert svd.config.svd_callback is svd_rayleigh_ritz_callback
+        svd.run_svd_analysis()
 
-        assert svd.config.svd_callback is svd_dense
+        # SVD Rayleigh-Ritz is not consistent with signs - manually iterate and check abs value
+        for i in range(5):
+            for j in range(4):
+                if (i, j) in [(1, 1), (2, 3), (3, 0), (4, 2)]:
+                    assert abs(svd.u[i, j]) == pytest.approx(1, abs=1e-6, rel=1e-6)
+                else:
+                    assert svd.u[i, j] == pytest.approx(0, abs=1e-6)
+
+        np.testing.assert_array_almost_equal(svd.s, np.array([0.1, 1, 5, 10]))
 
+        for i in range(5):
+            for j in range(4):
+                if (i, j) in [(1, 1), (2, 3), (3, 0), (4, 2)]:
+                    assert abs(svd.v[i, j]) == pytest.approx(1, abs=1e-6, rel=1e-6)
+                else:
+                    assert svd.v[i, j] == pytest.approx(0, abs=1e-6)
+
+    @pytest.mark.unit
+    def test_run_svd_analysis_dense(self, dummy_problem):
+        svd = SVDToolbox(dummy_problem, svd_callback=svd_dense)
         svd.run_svd_analysis()
 
         np.testing.assert_array_almost_equal(
@@ -2202,14 +2223,14 @@ def test_display_rank_of_equality_constraints(self, dummy_problem):
 
     @pytest.mark.unit
     def test_display_rank_of_equality_constraints(self, dummy_problem):
-        svd = SVDToolbox(dummy_problem, singular_value_tolerance=1)
+        svd = SVDToolbox(dummy_problem, singular_value_tolerance=0.9)
 
         stream = StringIO()
         svd.display_rank_of_equality_constraints(stream=stream)
 
         expected = """====================================================================================
 
-Number of Singular Values less than 1.0E+00 is 1
+Number of Singular Values less than 9.0E-01 is 1
 
 ====================================================================================
 """
@@ -2302,7 +2323,7 @@ def test_display_underdetermined_variables_and_constraints_specific(
 
     @pytest.mark.unit
     def test_display_underdetermined_variables_and_constraints(self, dummy_problem):
-        svd = SVDToolbox(dummy_problem, size_cutoff_in_singular_vector=1)
+        svd = SVDToolbox(dummy_problem, size_cutoff_in_singular_vector=1.1)
 
         stream = StringIO()
         svd.display_underdetermined_variables_and_constraints(stream=stream)