'Solve a linear system of equations with factor-structured
coefficients.
Parameters
rhs : array-like
A 2 dimensional array with the same number of rows as the
PSD matrix represented by the class instance.
Returns
C^{-1} * rhs, where C is the covariance matrix represented
by this class instance.
Notes
This function exploits the factor structure for efficiency.'

def solve(self, rhs):

    qval = ((- self.scales) / (1 + self.scales))
    dr = np.sqrt(self.diag)
    rhs = (rhs / dr[:, None])
    mat = (qval[:, None] * np.dot(self.factor.T, rhs))
    rhs = (rhs + np.dot(self.factor, mat))
    return (rhs / dr[:, None])
