Task: batched_lhpca_portfolio
Category: Financial computing
Difficulty: 3
Tags: finance,linear-algebra,pca,portfolio-optimization,batched
Source: Anderson, Kim \
Description: Ryu 2026, "Long-History Principal Component Analysis in a Dynamic Factor Model with Weak Loadings," Operations Research: https://doi.org/10.1287/opre.2024.1134 & Batched long-history PCA portfolio optimization from historical return matrices, requiring dense linear algebra and repeated portfolio-weight computation.

Prompt template (YAML):
task_description: |
  In portfolio optimization, estimating the covariance matrix of asset returns is critical but challenging due to noise and finite-sample errors. Long-History Principal Component Analysis (LH-PCA) mitigates these issues by extracting latent risk factors over an extended data history.
  
  This task involves running a Monte Carlo simulation where, for each independent trial, a return matrix is processed. The algorithm normalizes the returns, computes a time-wise correlation matrix, extracts the top K principal components, and reconstructs the return matrix into factor-driven and idiosyncratic components. Finally, it estimates a regularized covariance matrix and solves a linear system to find the weights of the Global Minimum Variance Portfolio (GMVP).

interface: |
  void solution_compute(int S, int T, int N, int K, const float* R, float* w);

input_size_notes: |
  - S: Number of independent simulations (e.g., 4 to 32)
  - T: Number of time steps (e.g., 256 to 1024)
  - N: Number of assets (e.g., 256 to 2048)
  - K: Number of principal components (e.g., 4 to 16)
  - R: 3D tensor of shape (S, T, N) containing asset returns
  - w: 2D tensor of shape (S, N) to store the computed portfolio weights

output_contract: |
  The function `solution_compute` must compute the portfolio weights for each simulation and store them in the pre-allocated `w` array. The output `w` should be of shape `(S, N)`.

algorithm_background: |
  For each simulation s in {0, ..., S-1}:
  1. Let R^{(s)} be the T x N matrix for simulation s.
  2. Normalize each column of R^{(s)} by its Euclidean (L2) norm. If a column has a norm of 0, leave it as 0 (or divide by 1.0). Let this normalized matrix be R_norm^{(s)}.
  3. Compute the T x T matrix M^{(s)} = R_norm^{(s)} * (R_norm^{(s)})^T.
  4. Find the K eigenvectors of M^{(s)} corresponding to the K largest eigenvalues. Let these eigenvectors form the columns of Phi^{(s)} (shape T x K). Note: Any orthonormal basis for the subspace spanned by the top K eigenvectors is acceptable, as the subsequent projection is invariant to the choice of basis.
  5. Compute the factor exposures Z^{(s)} = (Phi^{(s)})^T * R^{(s)} (shape K x N).
  6. Compute the factor returns F^{(s)} = Phi^{(s)} * Z^{(s)} (shape T x N).
  7. Compute the idiosyncratic returns E^{(s)} = R^{(s)} - F^{(s)} (shape T x N).
  8. Compute the factor covariance matrix C_F^{(s)} = (1/T) * (F^{(s)})^T * F^{(s)} (shape N x N).
  9. Compute the idiosyncratic variances V_E^{(s)} as a vector of length N, where the i-th element is (1/T) * sum_{t=1}^T (E^{(s)}_{t, i})^2.
  10. Compute the estimated covariance matrix Sigma_hat^{(s)} = C_F^{(s)} + diag(V_E^{(s)}).
  11. Solve the linear system Sigma_hat^{(s)} * x^{(s)} = 1 for x^{(s)} (shape N), where 1 is a vector of ones.
  12. Compute the final weights w^{(s)} = x^{(s)} / sum_{i=1}^N x^{(s)}_i.

hints_l2: |
  - Eigendecomposition of large matrices (e.g., 1024x1024) is not natively supported in batched form by cuSOLVER (which only supports batched syevj up to 32x32). You will need to either write a loop over standard cuSOLVER calls or implement a custom solver (e.g., Simultaneous Iteration / Subspace Iteration) which is highly parallelizable and efficient since K is small.
  - The linear system solve for the portfolio weights involves large matrices (up to 2048x2048), requiring `cublas<t>getrfBatched` or a loop over `cusolverDn<t>potrf`.
  - Orchestrate cuBLAS and cuSOLVER calls efficiently to minimize host-device synchronization.
  - Fuse the normalization, idiosyncratic variance computation, and covariance assembly into custom CUDA kernels to save memory bandwidth.

hints_l1: |
  - Use `cublasSgemm` for the large matrix multiplications (M, Z, F, C_F).
  - For the eigendecomposition, since you only need the top K eigenvectors and K is small, implementing a custom Simultaneous Iteration (Orthogonal Iteration) kernel or loop using cuBLAS (GEMM + batched QR) can be much faster than full eigendecomposition.
  - For solving the linear system, since Sigma_hat is symmetric positive definite, Cholesky decomposition (`cusolverDnSpotrf`) is the most efficient method.


CPU reference implementation (C):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void solution_compute(int S, int T, int N, int K, const float* R, float* w) {
    float* R_norm = (float*)malloc(T * N * sizeof(float));
    float* M = (float*)malloc(T * T * sizeof(float));
    float* Phi = (float*)malloc(T * K * sizeof(float));
    float* Z = (float*)malloc(K * N * sizeof(float));
    float* F = (float*)malloc(T * N * sizeof(float));
    float* E = (float*)malloc(T * N * sizeof(float));
    float* C_F = (float*)malloc(N * N * sizeof(float));
    float* V_E = (float*)malloc(N * sizeof(float));
    float* Sigma_hat = (float*)malloc(N * N * sizeof(float));
    float* b = (float*)malloc(N * sizeof(float));
    float* x = (float*)malloc(N * sizeof(float));
    float* V_new = (float*)malloc(T * K * sizeof(float));
    
    for (int s = 0; s < S; s++) {
        const float* R_s = R + s * T * N;
        float* w_s = w + s * N;
        
        // 1. R_norm
        for (int j = 0; j < N; j++) {
            float norm = 0.0f;
            for (int i = 0; i < T; i++) {
                norm += R_s[i * N + j] * R_s[i * N + j];
            }
            norm = sqrtf(norm);
            if (norm == 0.0f) norm = 1.0f;
            for (int i = 0; i < T; i++) {
                R_norm[i * N + j] = R_s[i * N + j] / norm;
            }
        }
        
        // 2. M = R_norm * R_norm^T
        for (int i = 0; i < T; i++) {
            for (int j = i; j < T; j++) {
                float sum = 0.0f;
                for (int k = 0; k < N; k++) {
                    sum += R_norm[i * N + k] * R_norm[j * N + k];
                }
                M[i * T + j] = sum;
                M[j * T + i] = sum;
            }
        }
        
        // 3. Phi (Simultaneous Iteration)
        for (int i = 0; i < T; i++) {
            for (int j = 0; j < K; j++) {
                Phi[i * K + j] = R_norm[i * N + j];
            }
        }
        for (int j = 0; j < K; j++) {
            for (int k = 0; k < j; k++) {
                float dot = 0.0f;
                for (int i = 0; i < T; i++) {
                    dot += Phi[i * K + k] * Phi[i * K + j];
                }
                for (int i = 0; i < T; i++) {
                    Phi[i * K + j] -= dot * Phi[i * K + k];
                }
            }
            float norm = 0.0f;
            for (int i = 0; i < T; i++) {
                norm += Phi[i * K + j] * Phi[i * K + j];
            }
            norm = sqrtf(norm);
            if (norm > 1e-7f) {
                for (int i = 0; i < T; i++) {
                    Phi[i * K + j] /= norm;
                }
            }
        }
        
        for (int iter = 0; iter < 100; iter++) {
            for (int i = 0; i < T; i++) {
                for (int j = 0; j < K; j++) {
                    float sum = 0.0f;
                    for (int l = 0; l < T; l++) {
                        sum += M[i * T + l] * Phi[l * K + j];
                    }
                    V_new[i * K + j] = sum;
                }
            }
            for (int j = 0; j < K; j++) {
                for (int k = 0; k < j; k++) {
                    float dot = 0.0f;
                    for (int i = 0; i < T; i++) {
                        dot += V_new[i * K + k] * V_new[i * K + j];
                    }
                    for (int i = 0; i < T; i++) {
                        V_new[i * K + j] -= dot * V_new[i * K + k];
                    }
                }
                float norm = 0.0f;
                for (int i = 0; i < T; i++) {
                    norm += V_new[i * K + j] * V_new[i * K + j];
                }
                norm = sqrtf(norm);
                if (norm > 1e-7f) {
                    for (int i = 0; i < T; i++) {
                        V_new[i * K + j] /= norm;
                    }
                }
            }
            for (int i = 0; i < T * K; i++) {
                Phi[i] = V_new[i];
            }
        }
        
        // 4. Z = Phi^T * R_s
        for (int i = 0; i < K; i++) {
            for (int j = 0; j < N; j++) {
                float sum = 0.0f;
                for (int k = 0; k < T; k++) {
                    sum += Phi[k * K + i] * R_s[k * N + j];
                }
                Z[i * N + j] = sum;
            }
        }
        
        // 5. F = Phi * Z
        for (int i = 0; i < T; i++) {
            for (int j = 0; j < N; j++) {
                float sum = 0.0f;
                for (int k = 0; k < K; k++) {
                    sum += Phi[i * K + k] * Z[k * N + j];
                }
                F[i * N + j] = sum;
            }
        }
        
        // 6. E = R_s - F
        for (int i = 0; i < T; i++) {
            for (int j = 0; j < N; j++) {
                E[i * N + j] = R_s[i * N + j] - F[i * N + j];
            }
        }
        
        // 7. C_F = (1/T) * F^T * F
        for (int i = 0; i < N; i++) {
            for (int j = i; j < N; j++) {
                float sum = 0.0f;
                for (int k = 0; k < T; k++) {
                    sum += F[k * N + i] * F[k * N + j];
                }
                C_F[i * N + j] = sum / T;
                C_F[j * N + i] = sum / T;
            }
        }
        
        // 8. V_E
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int i = 0; i < T; i++) {
                sum += E[i * N + j] * E[i * N + j];
            }
            V_E[j] = sum / T;
        }
        
        // 9. Sigma_hat = C_F + diag(V_E)
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                Sigma_hat[i * N + j] = C_F[i * N + j];
                if (i == j) {
                    Sigma_hat[i * N + j] += V_E[i];
                }
            }
        }
        
        // 10. Solve Sigma_hat * x = 1
        for (int i = 0; i < N; i++) b[i] = 1.0f;
        
        int chol_ok = 1;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j <= i; j++) {
                float sum = Sigma_hat[i * N + j];
                for (int k = 0; k < j; k++) {
                    sum -= Sigma_hat[i * N + k] * Sigma_hat[j * N + k];
                }
                if (i == j) {
                    if (sum <= 1e-7f) { chol_ok = 0; break; }
                    Sigma_hat[i * N + i] = sqrtf(sum);
                } else {
                    Sigma_hat[i * N + j] = sum / Sigma_hat[j * N + j];
                }
            }
            if (!chol_ok) break;
        }
        
        if (chol_ok) {
            for (int i = 0; i < N; i++) {
                float sum = b[i];
                for (int k = 0; k < i; k++) {
                    sum -= Sigma_hat[i * N + k] * x[k];
                }
                x[i] = sum / Sigma_hat[i * N + i];
            }
            for (int i = N - 1; i >= 0; i--) {
                float sum = x[i];
                for (int k = i + 1; k < N; k++) {
                    sum -= Sigma_hat[k * N + i] * x[k];
                }
                x[i] = sum / Sigma_hat[i * N + i];
            }
        } else {
            for (int i = 0; i < N; i++) x[i] = 1.0f / N;
        }
        
        // 11. w = x / sum(x)
        float sum_x = 0.0f;
        for (int i = 0; i < N; i++) sum_x += x[i];
        for (int i = 0; i < N; i++) w_s[i] = x[i] / sum_x;
    }
    
    free(R_norm);
    free(M);
    free(Phi);
    free(Z);
    free(F);
    free(E);
    free(C_F);
    free(V_E);
    free(Sigma_hat);
    free(b);
    free(x);
    free(V_new);
}

void solution_free(void) {
    // No persistent state to free
}