Math.js regression performance problem, code

output.txt

_xMatrixTransp size =  2,4084
_xMatrix size =  4084
_normalTargetData size =  4084

regression-thread-code.js

function numericJs() {
  var xMatrix = numeric.transpose(xMatrixTransp);
  // see https://en.wikipedia.org/wiki/Ordinary_least_squares#Estimation
  // Compute beta = (X^T X)^{-1} X^T y 
  var dotProduct = numeric.dot(xMatrixTransp, xMatrix);
  var inverse = numeric.inv(dotProduct);
  dotProduct = null;
  var multi2 = numeric.dot(inverse, xMatrixTransp);
  var betas = numeric.dot(multi2, normalTargetData);
  multi2 = null;

  // console.log("numericjs betas", betas);

  var n = _.size(xMatrix),
  k = _.size(xMatrix[0]),
  beta = betas[1];

  // console.log("n, k", n, k);

  // get confidence interval and p-value
  var ciAndPvalue = regressionUtils.getCIAndPvalueNumericjs(inverse, xMatrix, xMatrixTransp, [normalTargetData], n, k, beta);

  return {
    result: { success: true },
    ci: ciAndPvalue.ci,
    pvalue: ciAndPvalue.pvalue,
    betas: betas
  };

}

function mathJs() {
  // see https://en.wikipedia.org/wiki/Ordinary_least_squares#Estimation
  // beta = (X^T X)^{-1} X^T y 
  var _xMatrixTransp = math.matrix(xMatrixTransp);
  var _xMatrix = math.transpose(xMatrixTransp);
  var _normalTargetData = math.matrix(normalTargetData);

  console.log("_xMatrixTransp size = ", String(_xMatrixTransp.size()));
  console.log("_xMatrix size = ", _.size(_xMatrix));
  console.log("_normalTargetData size = ", String(_normalTargetData.size()));

  // Compute beta = (X^T X)^{-1} X^T y 
  var dotProduct = math.multiply(_xMatrixTransp, _xMatrix);
  var inverse = math.inv(dotProduct);
  var multi2 = math.multiply(inverse, _xMatrixTransp);
  var betas = math.multiply(multi2, normalTargetData);

  // console.log("mathjs betas=", String(betas));

  var n = math.size(_xMatrix)[0];
  var k = math.size(_xMatrix)[1];
  var beta = betas.subset(math.index(1));

  var _yMatrixTransp = math.matrix([normalTargetData]);

  // get confidence interval and p-value
  var ciAndPvalue = regressionUtils.getCIAndPvalueMathjs(inverse, _xMatrix, _xMatrixTransp, _yMatrixTransp, n, k, beta);

  return {
    result: { success: true },
    ci: ciAndPvalue.ci,
    pvalue: ciAndPvalue.pvalue,
    betas: betas.valueOf()
  };

}
root.getCIAndPvalueMathjs = function(dotInverse, xMatrix, xMatrixTransposed, yMatrixTransposed, n, k, beta) {
  var VARIABLE_INDEX = 1;

  // H = X (X^T X)^{-1} X^T
  var hMatrix = math.chain(xMatrix)
  .multiply(dotInverse)
  .multiply(xMatrixTransposed)
  .done();

  xMatrix = null;

  var hMatrixSize = math.size(hMatrix).subset(math.index(0));
  var identity = math.eye(hMatrixSize, 'sparse');
  var subtracted = math.subtract(identity, hMatrix);
  var yMatrix = math.transpose(yMatrixTransposed);

  var degrees = n-(k+1);

  var sigma = math.chain(yMatrixTransposed)
  .multiply(subtracted)
  .multiply(yMatrix)
  .divide(degrees)
  .done();

  yMatrix = null;
  subtracted = null;
  identity = null;
  yMatrixTransposed = null;

  var cMatrix = math.multiply(sigma, dotInverse);

  dotInverse = null;

  var alpha = 0.05 / k;

  var sqrt = Math.sqrt( cMatrix.subset(math.index(VARIABLE_INDEX, VARIABLE_INDEX)) );
  cMatrix = null;
  var ci = statDist.tdistr(degrees, alpha/2) * sqrt;

  // See http://reliawiki.org/index.php/Multiple_Linear_Regression_Analysis 
  // -> p value = 2 * (1-P(T <= |t0|)
  var t = beta / sqrt;
  sqrt = null;
  var pvalue = statDist.tprob(degrees, t);

  return {
    ci: [ beta - ci, beta + ci ],
    pvalue: pvalue
  };
};

root.getCIAndPvalueNumericjs = function(dotInverse, xMatrix, xMatrixTransposed, yMatrixTransposed, n, k, beta) {
    var VARIABLE_INDEX = 1;

    var hMatrix = numeric.dot( numeric.dot(xMatrix, dotInverse), xMatrixTransposed );

    var _identity = numeric.identity( _.size(hMatrix) );

    var _subtracted = numeric.sub(_identity, hMatrix);

    var yMatrix = numeric.transpose(yMatrixTransposed);

    var sigma = numeric.dot( numeric.dot( yMatrixTransposed, _subtracted ), yMatrix )[0][0] / (n-(k+1));
    var cMatrix = numeric.mul(sigma, dotInverse);

    var alpha = 0.05 / k;

    var _sqrt = Math.sqrt( cMatrix[VARIABLE_INDEX][VARIABLE_INDEX] ); //numeric.getDiag(cMatrix)[1] );
    var degrees = n-(k+1);
    var ci = statDist.tdistr(degrees, alpha/2) * _sqrt;

    // See http://reliawiki.org/index.php/Multiple_Linear_Regression_Analysis 
    // -> p value = 2 * (1-P(T <= |t0|)
    var t = beta / _sqrt;
    var pvalue = statDist.tprob(degrees, t);

    return {
      ci: [ beta - ci, beta + ci ],
      pvalue: pvalue
  };
};