robust_solve_XtX function

Robust solve for X'X beta = X'Y using QR decomposition