Refactor for Performance and Clarity

A quant on the desk wrote this code to build a pairwise correlation matrix across a basket of tickers' daily return series, so it can feed into a portfolio risk model. It's correct on the small baskets used in backtests, but on the full universe of a few hundred tickers with years of daily returns it's become too slow to run intraday.

double computeCorrelation(std::vector<double> a, std::vector<double> b) {
    double sumA = 0, sumB = 0;
    for (int i = 0; i < a.size(); i++) {
        sumA += a[i];
        sumB += b[i];
    }
    double meanA = sumA / a.size();
    double meanB = sumB / b.size();

    double cov = 0, varA = 0, varB = 0;
    for (int i = 0; i < a.size(); i++) {
        cov += (a[i] - meanA) * (b[i] - meanB);
        varA += (a[i] - meanA) * (a[i] - meanA);
        varB += (b[i] - meanB) * (b[i] - meanB);
    }

    return cov / sqrt(varA * varB);
}

std::vector<std::vector<double>> buildCorrelationMatrix(std::vector<std::vector<double>> returns) {
    std::vector<std::vector<double>> matrix;
    for (int i = 0; i < returns.size(); i++) {
        std::vector<double> row;
        for (int j = 0; j < returns.size(); j++) {
            row.push_back(computeCorrelation(returns[i], returns[j]));
        }
        matrix.push_back(row);
    }
    return matrix;
}

Refactor this for performance while keeping it correct. Explain exactly what makes the current version slow and what could go wrong with it beyond just speed, before proposing your changes.

Implement correlation_matrix(returns), where returns is a nonempty list of equal-length numeric series. Return a square list of floating-point correlations. Treat a zero-variance series as having correlation 0.0 with every series, including itself.