Z - MxQ
    mu - NxQ
    S - NxQ
    gamma - NxQ

def _psi2computations(dL_dpsi2, variance, Z, mu, S, gamma):
    """
    Z - MxQ
    mu - NxQ
    S - NxQ
    gamma - NxQ
    """
    # here are the "statistics" for psi1 and psi2
    # Produced intermediate results:
    # _psi2_dvariance      Q
    # _psi2_dZ             MxQ
    # _psi2_dgamma         NxQ
    # _psi2_dmu            NxQ
    # _psi2_dS             NxQ
    
    mu2 = np.square(mu)
    gamma2 = np.square(gamma)
    variance2 = np.square(variance)
    mu2S = mu2+S # NxQ
    gvm = gamma*mu*variance
    common_sum = gvm.dot(Z.T)
    Z_expect = (np.dot(dL_dpsi2,Z)*Z).sum(axis=0)
    Z_expect_var2 = Z_expect*variance2
    dL_dpsi2T = dL_dpsi2+dL_dpsi2.T
    common_expect = common_sum.dot(dL_dpsi2T).dot(Z)
    Z2_expect = common_sum.dot(dL_dpsi2T)
    Z1_expect = dL_dpsi2T.dot(Z)
    
    dL_dvar = variance*Z_expect*2.*(gamma*mu2S-gamma2*mu2).sum(axis=0)+(common_expect*gamma*mu).sum(axis=0)
        
    dL_dgamma = Z_expect_var2*(mu2S-2.*gamma*mu2)+common_expect*mu*variance
                
    dL_dmu = Z_expect_var2*mu*2.*(gamma-gamma2) + common_expect*gamma*variance

    dL_dS = gamma*Z_expect_var2
    
    dL_dZ = (gamma*(mu2S-gamma*mu2)).sum(axis=0)*variance2*Z1_expect+ Z2_expect.T.dot(gamma*mu)*variance

    return dL_dvar, dL_dgamma, dL_dmu, dL_dS, dL_dZ