Skip to content

Commit

Permalink
Lambda is now vectorized when provided
Browse files Browse the repository at this point in the history
Looting wire transfer xecutive SERT offensive information warfare
counter terrorism SEAL Team 6 sweep Typhoon Ruby Ridge Planet-1
remailers H1N1 Food Poisoning Pine Gap
  • Loading branch information
benjamin-james committed Mar 14, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
1 parent b1f8258 commit 062ee11
Showing 1 changed file with 11 additions and 1 deletion.
12 changes: 11 additions & 1 deletion src/se.hpp
Original file line number Diff line number Diff line change
@@ -49,6 +49,16 @@ Eigen::ArrayXd robust_se_X(const Eigen::MatrixBase<TX> &Xmat,
return tval * var.cwiseMax(epsilon).array().rsqrt();
}

template<typename TY, typename TL>
Eigen::VectorXd ols_beta1(Eigen::VectorXd X, double X_sqnorm,
const Eigen::MatrixBase<TY> &Y,
const Eigen::ArrayBase<TL> &lambda)
{
Eigen::ArrayXd numer = (X.transpose() * Y).eval().array();
Eigen::ArrayXd denom = (lambda + X_sqnorm).eval();
Eigen::ArrayXd quot = numer / denom;
return quot.matrix();
}
template<typename TX, typename TY, typename TL>
Eigen::ArrayXd robust_se_L(const Eigen::MatrixBase<TX> &Xmat,
const Eigen::MatrixBase<TY> &Y, /* V\Sigma?? */
@@ -57,8 +67,8 @@ Eigen::ArrayXd robust_se_L(const Eigen::MatrixBase<TX> &Xmat,
{
Eigen::VectorXd X = Xmat.col(0).eval();
double X_sqnorm = std::max(X.squaredNorm(), epsilon);
Eigen::VectorXd beta = ols_beta1(X, X_sqnorm, Y, lambda);
// X pseudoinverse is X' / squared_norm
Eigen::VectorXd beta = ((X.transpose() * Y).array() / (X_sqnorm + lambda)).matrix().eval();
Eigen::VectorXd var = ((X.transpose() / X_sqnorm).cwiseAbs2() * (Y - X * beta.transpose()).cwiseAbs2()).eval();
Eigen::ArrayXd tval = beta.array();
return tval * var.cwiseMax(epsilon).array().rsqrt();

0 comments on commit 062ee11

Please sign in to comment.