import numpy as np

def compute_symmetric_matrix_sqrt(M: np.ndarray) -> np.ndarray:
    """
    Computes the square root of a symmetric PSD matrix using eigen-decomposition.

    Parameters:
    - M (np.ndarray): A symmetric positive semidefinite matrix.

    Returns:
    - np.ndarray: The matrix square root of H_inv.
    """
    eigvals, eigvecs = np.linalg.eigh(M)
    sqrt_eigvals = np.sqrt(np.clip(eigvals, a_min=0, a_max=None))  # avoid tiny negatives
    sqrt_M = eigvecs @ np.diag(sqrt_eigvals) @ eigvecs.T
    return sqrt_M