'Solves matrix equations of the form `covmat * soln = rhs` and
returns the values of `soln`, where `covmat` is the covariance
matrix represented by this class.
Parameters
expval: array-like
The expected value of endog for each observed value in the
group.
index: integer
The group index.
stdev : array-like
The standard deviation of endog for each observation in
the group.
rhs : list/tuple of array-like
A set of right-hand sides; each defines a matrix equation
to be solved.
Returns
soln : list/tuple of array-like
The solutions to the matrix equations.
Notes
Returns None if the solver fails.
Some dependence structures do not use `expval` and/or `index`
to determine the correlation matrix.  Some families
(e.g. binomial) do not use the `stdev` parameter when forming
the covariance matrix.
If the covariance matrix is singular or not SPD, it is
projected to the nearest such matrix.  These projection events
are recorded in the fit_history member of the GEE model.
Systems of linear equations with the covariance matrix as the
left hand side (LHS) are solved for different right hand sides
(RHS); the LHS is only factorized once to save time.
This is a default implementation, it can be reimplemented in
subclasses to optimize the linear algebra according to the
struture of the covariance matrix.'

def covariance_matrix_solve(self, expval, index, stdev, rhs):

    (vmat, is_cor) = self.covariance_matrix(expval, index)
    if is_cor:
        vmat *= np.outer(stdev, stdev)
    threshold = 0.01
    success = False
    cov_adjust = 0
    for itr in range(20):
        try:
            vco = spl.cho_factor(vmat)
            success = True
            break
        except np.linalg.LinAlgError:
            vmat = cov_nearest(vmat, method=self.cov_nearest_method, threshold=threshold)
            threshold *= 2
            cov_adjust += 1
    self.cov_adjust.append(cov_adjust)
    if (not success):
        warnings.warn('Unable    to    condition    covariance    matrix    to    an    SPD    matrix    using    cov_nearest', ConvergenceWarning)
        vmat = np.diag(np.diag(vmat))
        vco = spl.cho_factor(vmat)
    soln = [spl.cho_solve(vco, x) for x in rhs]
    return soln
