def _smw_solver(s, A, AtA, BI, di):

    # Use SMW identity
    qmat = AtA / s
    m = BI.shape[0]
    qmat[0:m, 0:m] += BI
    ix = np.arange(m, A.shape[1])
    qmat[ix, ix] += di
    if sparse.issparse(A):
        qi = sparse.linalg.inv(qmat)
        qmati = A.dot(qi.T).T
    else:
        qmati = np.linalg.solve(qmat, A.T)

    def solver(rhs):
        if sparse.issparse(A):
            ql = qmati.dot(rhs)
            ql = A.dot(ql)
        else:
            ql = np.dot(qmati, rhs)
            ql = np.dot(A, ql)
        rslt = rhs / s - ql / s**2
        if sparse.issparse(rslt):
            rslt = np.asarray(rslt.todense())
        return rslt

    return solver
Solves the system (s*I + A*B*A') * x = rhs for an arbitrary rhs.
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/regression/mixed_linear_model.py/_smw_solver