def METHOD_NAME(
    ideal_matrix: npt.NDArray[np.complex64],
    basis_matrices: List[npt.NDArray[np.complex64]],
    tol: float = 1.0e-8,
    initial_guess: Optional[npt.NDArray[np.float64]] = None,
) -> npt.NDArray[np.float64]:
    r"""
    Returns the list of real coefficients :math:`[x_0, x_1, \dots]`,
    which minimizes :math:`\sum_j |x_j|` with the contraint that
    the following representation of the input ``ideal_matrix`` holds:
    .. math::
        \text{ideal_matrix} = x_0 A_0 + x_1 A_1 + \cdots
    where :math:`\{A_j\}` are the basis matrices, i.e., the elements of
    the input ``basis_matrices``.
    This function can be used to compute the optimal representation
    of an ideal superoperator (or Choi state) as a linear
    combination of real noisy superoperators (or Choi states).
    Args:
        ideal_matrix: The ideal matrix to represent.
        basis_matrices: The list of basis matrices.
        tol: The error tolerance for each matrix element
            of the represented matrix.
        initial_guess: Optional initial guess for the coefficients
            :math:`[x_0, x_1, \dots]`.
    Returns:
        The list of optimal coefficients :math:`[x_0, x_1, \dots]`.
    """
    # Map complex matrices to extended real matrices
    ideal_matrix_real = np.hstack(
        (np.real(ideal_matrix), np.imag(ideal_matrix))
    )
    basis_matrices_real = [
        np.hstack((np.real(mat), np.imag(mat))) for mat in basis_matrices
    ]
    # Express the representation constraint written in the docstring in the
    # form of a matrix multiplication applied to the x vector: A @ x == b.
    matrix_a = np.array(
        [
            matrix_to_vector(mat)  # type: ignore[arg-type]
            for mat in basis_matrices_real
        ]
    ).T
    array_b = matrix_to_vector(ideal_matrix_real)  # type: ignore[arg-type]
    constraint = LinearConstraint(matrix_a, lb=array_b - tol, ub=array_b + tol)
    def one_norm(x: npt.NDArray[np.complex64]) -> float:
        return cast(float, np.linalg.norm(x, 1))
    if initial_guess is None:
        initial_guess = np.zeros(len(basis_matrices))
    result = minimize(one_norm, x0=initial_guess, constraints=constraint)
    if not result.success:
        raise RuntimeError("The search for an optimal representation failed.")
    return result.x