diff --git a/DESCRIPTION b/DESCRIPTION index ece5dd9..eeb89ae 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: gputools -Version: 2.1.0 +Version: 2.1.1 Title: GPU Drop-Ins for Statistics on Wide Data Authors@R: c(person("Josh", "Buckner", email = "nullsatz@gmail.com", role = c("aut", "cre")), diff --git a/NEWS.md b/NEWS.md index 44aa18a..7d73779 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,20 @@ +# gputools 2.1.1 + +- `gpuGranger()` reported F = 0, p = 1, "no causality", for a pair it + could not resolve (collinear lags, a constant series, a cause identical + to the effect), computed in single precision by default on uncentred + data (F statistics of long random walks off by more than 10 percent), + and built one regression per pair with memory growing as the square of + the number of series. It now computes every pair in double by the + Frisch-Waugh route, with each effect's own past projected out once, and + the cause's test taken from a few matrix products, agreeing with `lm()` + to about 1e-8; a pair it cannot resolve in double is `NA` with a warning + (series integrated twice or more should be differenced first), a series + against itself or a constant series `NA` without one. `lag` is checked + (1 to 16). The `precision` argument is accepted but the tests always + compute in double. The method and its limits follow the Python package + widestats. + # gputools 2.1.0 - New `gpuCorPerm()`: permutation p-values for every correlation between diff --git a/R/gpuGranger.R b/R/gpuGranger.R index 71227d9..bb2b916 100644 --- a/R/gpuGranger.R +++ b/R/gpuGranger.R @@ -1,21 +1,31 @@ gpuGranger <- function(x, y=NULL, lag, precision = getOption("gputools.precision", "single")) { + # the argument is accepted for the common interface; the Granger tests + # always compute in double (see the help page) precision <- .precision(precision) x <- as.matrix(x) rows <- nrow(x) colsx <- ncol(x) - if(rows - lag <= 2 * lag + 1) { + if (ncol(x) < 1) stop("x needs at least one column") + if (length(lag) != 1 || !is.numeric(lag) || is.na(lag) || + lag != round(lag) || lag < 1) { + stop("lag must be a positive whole number") + } + if (lag > 16) { + stop("lag must be at most 16") + } + lag <- as.integer(lag) + if (rows - lag <= 2 * lag + 1) { stop("time sequence too short for lag: use longer sequences or smaller lag") } if(!is.null(y)) { y <- as.matrix(y) + if (nrow(y) != rows) stop("x and y have different numbers of rows") } - lag <- as.integer(lag) - if(is.null(y)) { colsy <- colsx cRetVal <- .Call(C_gpuGranger, x, as.integer(rows), as.integer(colsx), @@ -26,6 +36,14 @@ gpuGranger <- function(x, y=NULL, lag, as.integer(colsx), y, as.integer(colsy), lag, precision) } + if (cRetVal$unresolved > 0) { + warning(sprintf(paste( + "%d pair(s) are too ill-conditioned to resolve in double", + "precision: a cause almost exactly a linear function of the", + "effect's own past, or nearly collinear lags (a series integrated", + "twice or more: difference it; or an exactly periodic one); their", + "tests are NA"), cRetVal$unresolved), call. = FALSE) + } # [i, j] tests whether x[, i] Granger-causes x[, j] (or y[, j]); the # dimnames say so when printed names <- function(m, prefix) @@ -35,10 +53,9 @@ gpuGranger <- function(x, y=NULL, lag, dn <- list(cause = causes, effect = effects) fStats <- matrix(cRetVal$fStats, colsx, colsy, dimnames = dn) pValues <- matrix(cRetVal$pValues, colsx, colsy, dimnames = dn) - if (is.null(y)) { - # a series does not Granger-cause itself: no test, so no value - diag(fStats) <- NA - diag(pValues) <- NA - } + # no test (a series against itself, a constant series) or no answer + # (unresolvable) is NA, not NaN, as R reports a missing statistic + fStats[is.nan(fStats)] <- NA + pValues[is.nan(pValues)] <- NA return(list(fStatistics = fStats, pValues = pValues)) } diff --git a/inst/cuda/granger.cu b/inst/cuda/granger.cu index df8cad8..6584b74 100644 --- a/inst/cuda/granger.cu +++ b/inst/cuda/granger.cu @@ -1,210 +1,276 @@ -#define max(a, b) ((a > b)?a:b) +// Granger causality F tests for every ordered pair of series, by the +// Frisch-Waugh route the host code in src/granger.cpp drives: +// +// - each effect's restricted design D = [1, own lags] (T x m, m = lag + 1) +// is orthonormalised to Q by Cholesky-QR2 (the kernels here factor the +// small Gram matrices and apply R^{-1}; the Gram products and the big +// products below are cuBLAS calls from the host); +// - with e the effect's residual and L every cause's lags (T x p*lag), +// the gain RSS_r - RSS_u of a pair is g' G^{-1} g for the cause's lags +// with the effect's own past projected out, M = L_i - Q Q' L_i: G = M'M, +// g = M'e; +// - the fast form computes those from the small products C = Q'L, L'e and +// Q'e as G = L'L - C'C and g = L'e - C'(Q'e), and a per-pair score says +// whether that cancellation can be trusted; the pairs it cannot be +// trusted on are recomputed from M formed explicitly. +// +// Everything is in double (the host always launches the double build); +// lag is at most GRANGER_MAXLAG so the per-pair systems fit in registers. +// All matrices are column-major, as cuBLAS has them. -#define THREADSPERDIM 16 +#define GRANGER_MAXLAG 16 +#define GRANGER_MAXM (GRANGER_MAXLAG + 1) -#define FALSE 0 -#define TRUE !FALSE - -// mX has order rows x cols -// vectY has length rows -__global__ void getRestricted(int countx, int county, int rows, int cols, - real * mX, int mXdim, real * vY, int vYdim, real * mQ, int mQdim, - real * mR, int mRdim, real * vectB, int vectBdim) { - - int - m = blockIdx.x * THREADSPERDIM + threadIdx.x, n, - i, j, k; - real - sum, invnorm, - * X, * Y, * Q, * R, * B, - * coli, * colj, - * colQ, * colX; - - if(m >= county) return; - // Any pair (m, n) holds column m's own lags; the single-matrix form - // leaves the diagonal pair (m, m) unfilled, so avoid n == m there, and the - // x/y form with one x column only has n == 0. - if(m == 1 || countx == 1) n = 0; - else n = 1; - - X = mX + (m * mXdim); - // initialize the intercepts - for(i = 0; i < rows; i++) - X[i] = real(1); - - Y = vY + (m * countx + n) * vYdim; - B = vectB + m * vectBdim; - Q = mQ + m * mQdim; - R = mR + m * mRdim; - - // initialize Q with X ... - for(i = 0; i < rows; i++) { - for(j = 0; j < cols; j++) - Q[i+j*rows] = X[i+j*rows]; +// Cholesky of the symmetric n x n matrix a (column-major, upper triangle +// read) into the upper factor r (a = r' r). Returns 0, or k > 0 if pivot k +// was not positive (or not finite). +__device__ int grangerCholesky(int n, const real * a, int lda, real * r) +{ + for (int j = 0; j < n; j++) { + for (int i = 0; i <= j; i++) { + real s = a[i + j * lda]; + for (int k = 0; k < i; k++) s -= r[k + i * n] * r[k + j * n]; + if (i == j) { + if (!(s > real(0)) || !isfinite(s)) return j + 1; + r[j + j * n] = sqrt(s); + } else { + r[i + j * n] = s / r[i + i * n]; + } + } + for (int i = j + 1; i < n; i++) r[i + j * n] = real(0); } + return 0; +} - // gramm-schmidt process to find Q - for(j = 0; j < cols; j++) { - colj = Q+rows*j; - for(i = 0; i < j; i++) { - coli = Q+rows*i; - sum = real(0); - for(k = 0; k < rows; k++) - sum += coli[k] * colj[k]; - for(k = 0; k < rows; k++) - colj[k] -= sum * coli[k]; +// Inverse of the upper triangular n x n matrix r into ri. +__device__ void grangerInvertUpper(int n, const real * r, real * ri) +{ + for (int j = 0; j < n; j++) { + for (int i = n - 1; i >= 0; i--) { + if (i > j) { ri[i + j * n] = real(0); continue; } + real s = (i == j) ? real(1) : real(0); + for (int k = i + 1; k <= j; k++) s -= r[i + k * n] * ri[k + j * n]; + ri[i + j * n] = s / r[i + i * n]; } - sum = real(0); - for(i = 0; i < rows; i++) - sum += colj[i] * colj[i]; - invnorm = real(1) / sqrt(sum); - for(i = 0; i < rows; i++) - colj[i] *= invnorm; } - for(i = cols-1; i > -1; i--) { - colQ = Q+i*rows; - // matmult Q * X -> R - for(j = 0; j < cols; j++) { - colX = X+j*rows; - sum = real(0); - for(k = 0; k < rows; k++) - sum += colQ[k] * colX[k]; - R[i+j*cols] = sum; - } - sum = real(0); - // compute the vector Q^t * Y -> B - for(j = 0; j < rows; j++) - sum += colQ[j] * Y[j]; - // back substitution to find the x for Rx = B - for(j = cols-1; j > i; j--) - sum -= R[i+j*cols] * B[j]; +} - B[i] = sum / R[i+i*cols]; - } +// ||r||_F ||r^{-1}||_F: at least the 2-norm condition number of r and at +// most n times it. +__device__ real grangerCondUpper(int n, const real * r, const real * ri) +{ + real a = real(0), b = real(0); + for (int j = 0; j < n; j++) + for (int i = 0; i <= j; i++) { + a += r[i + j * n] * r[i + j * n]; + b += ri[i + j * n] * ri[i + j * n]; + } + return sqrt(a) * sqrt(b); } -// mX has order rows x cols -// vectY has length rows -__global__ void getUnrestricted(int countx, int county, int rows, int cols, - real * mX, int mXdim, real * vY, int vYdim, real * mQ, int mQdim, - real * mR, int mRdim, real * vectB, int vectBdim) { +// The restricted design of each of b effects, side by side: D (T x b*m), +// effect jj's columns jj*m .. jj*m + m - 1 holding the intercept at unit +// length and its lag columns from ly (T x q*lag; effect j0 + jj's are +// (j0 + jj)*lag .. + lag - 1). +__global__ void grangerDesign(int T, int lag, int j0, int b, const real * ly, + real * d) +{ + int m = lag + 1; + size_t e = blockIdx.x * (size_t) blockDim.x + threadIdx.x; + if (e >= (size_t) T * b * m) return; + int t = e % T, col = e / T, jj = col / m, a = col % m; + d[e] = a == 0 ? rsqrt((real) T) + : ly[t + ((size_t) (j0 + jj) * lag + a - 1) * T]; +} - int - n = blockIdx.x * THREADSPERDIM + threadIdx.x, - m = blockIdx.y * THREADSPERDIM + threadIdx.y, - i, j, k; - real - sum, invnorm, - * X, * Y, * Q, * R, * B, - * coli, * colj, - * colQ, * colX; - if((m >= county) || (n >= countx)) return; +// From each effect's m x m Gram matrix g (= D'D, b of them): the inverse +// of its upper Cholesky factor into rinv, the factor's condition estimate +// into cond (on the first pass only, when writeCond is set), and bad set +// where the factorisation failed or the estimate is not below limit. +__global__ void grangerFactor(int b, int m, const real * g, real limit, + int writeCond, real * rinv, real * cond, + int * bad) +{ + int jj = blockIdx.x * blockDim.x + threadIdx.x; + if (jj >= b) return; + real r[GRANGER_MAXM * GRANGER_MAXM], ri[GRANGER_MAXM * GRANGER_MAXM]; + int info = grangerCholesky(m, g + (size_t) jj * m * m, m, r); + real c; + if (info) { + c = REAL_INF; + for (int k = 0; k < m * m; k++) ri[k] = REAL_NAN; + } else { + grangerInvertUpper(m, r, ri); + c = grangerCondUpper(m, r, ri); + } + if (writeCond) cond[jj] = c; + if (info || !(c < limit)) bad[jj] = 1; + for (int k = 0; k < m * m; k++) rinv[(size_t) jj * m * m + k] = ri[k]; +} - X = mX + (m * countx + n) * mXdim; - // initialize the intercepts - for(i = 0; i < rows; i++) - X[i] = real(1); +// Q = Q R^{-1} in place, each thread one row of one effect's T x m block. +__global__ void grangerApplyRinv(int T, int m, int b, const real * rinv, + real * q) +{ + size_t e = blockIdx.x * (size_t) blockDim.x + threadIdx.x; + if (e >= (size_t) T * b) return; + int t = e % T, jj = e / T; + real * row = q + (size_t) jj * m * T + t; + const real * ri = rinv + (size_t) jj * m * m; + real v[GRANGER_MAXM]; + for (int a = 0; a < m; a++) v[a] = row[(size_t) a * T]; + for (int c = 0; c < m; c++) { + real s = real(0); + for (int a = 0; a <= c; a++) s += v[a] * ri[a + c * m]; + row[(size_t) c * T] = s; + } +} - Y = vY + (m*countx+n) * vYdim; - B = vectB + (m*countx+n) * vectBdim; - Q = mQ + (m*countx+n) * mQdim; - R = mR + (m*countx+n) * mRdim; +// Each effect's residual sum of squares from its residual column of e. +__global__ void grangerRss(int T, int b, const real * e, real * rss) +{ + int jj = blockIdx.x * blockDim.x + threadIdx.x; + if (jj >= b) return; + const real * col = e + (size_t) jj * T; + real s = real(0); + for (int t = 0; t < T; t++) s += col[t] * col[t]; + rss[jj] = s; +} - // initialize Q with X ... - for(i = 0; i < rows; i++) { - for(j = 0; j < cols; j++) - Q[i+j*rows] = X[i+j*rows]; - } +// Each cause's lag columns' squared lengths: l2 (p*lag). +__global__ void grangerColumnSq(int T, int ncol, const real * l, real * l2) +{ + int c = blockIdx.x * blockDim.x + threadIdx.x; + if (c >= ncol) return; + const real * col = l + (size_t) c * T; + real s = real(0); + for (int t = 0; t < T; t++) s += col[t] * col[t]; + l2[c] = s; +} - // gramm-schmidt process to find Q - for(j = 0; j < cols; j++) { - colj = Q+rows*j; - for(i = 0; i < j; i++) { - coli = Q+rows*i; - sum = real(0); - for(k = 0; k < rows; k++) - sum += coli[k] * colj[k]; - for(k = 0; k < rows; k++) - colj[k] -= sum * coli[k]; - } - sum = real(0); - for(i = 0; i < rows; i++) - sum += colj[i] * colj[i]; - invnorm = real(1) / sqrt(sum); - for(i = 0; i < rows; i++) - colj[i] *= invnorm; +// Solves G x = g for a pair (G lag x lag, symmetric positive definite if +// the pair is resolvable) and returns the gain g'x, the Frobenius +// estimate of G's condition number (infinite if the factorisation fails, +// NaN if G is not finite) and the smallest diagonal fraction left. +__device__ real grangerSolvePair(int lag, const real * G, const real * g, + const real * l2i, real * cond, real * left) +{ + real r[GRANGER_MAXLAG * GRANGER_MAXLAG], ri[GRANGER_MAXLAG * GRANGER_MAXLAG]; + real lft = REAL_INF; + int finite = 1; + for (int k = 0; k < lag * lag; k++) if (!isfinite(G[k])) finite = 0; + for (int k = 0; k < lag; k++) { + real f = G[k + k * lag] / l2i[k]; + if (isnan(f)) { lft = f; break; } // a NaN propagates + if (f < lft) lft = f; } - for(i = cols-1; i > -1; i--) { - colQ = Q+i*rows; - // matmult Q * X -> R - for(j = 0; j < cols; j++) { - colX = X+j*rows; - sum = real(0); - for(k = 0; k < rows; k++) - sum += colQ[k] * colX[k]; - R[i+j*cols] = sum; + *left = lft; + if (!finite) { *cond = REAL_NAN; return REAL_NAN; } + if (grangerCholesky(lag, G, lag, r)) { *cond = REAL_INF; return REAL_NAN; } + grangerInvertUpper(lag, r, ri); + // ||G||_F ||G^{-1}||_F, G^{-1} = R^{-1} R^{-T} + real a = real(0), bsum = real(0); + for (int k = 0; k < lag * lag; k++) a += G[k] * G[k]; + for (int i = 0; i < lag; i++) + for (int j = 0; j < lag; j++) { + real s = real(0); + for (int k = (i > j ? i : j); k < lag; k++) s += ri[i + k * lag] * ri[j + k * lag]; + bsum += s * s; } - sum = real(0); - // compute the vector Q^t * Y -> B - for(j = 0; j < rows; j++) - sum += colQ[j] * Y[j]; - // back substitution to find the x for Rx = B - for(j = cols-1; j > i; j--) - sum -= R[i+j*cols] * B[j]; - - B[i] = sum / R[i+i*cols]; + *cond = sqrt(a) * sqrt(bsum); + // x = R^{-1} R^{-T} g, gain = g' x + real y[GRANGER_MAXLAG], gain = real(0); + for (int i = 0; i < lag; i++) { // R' y = g (forward) + real s = g[i]; + for (int k = 0; k < i; k++) s -= r[k + i * lag] * y[k]; + y[i] = s / r[i + i * lag]; } + for (int i = 0; i < lag; i++) gain += y[i] * y[i]; // g' G^{-1} g = y' y + return gain; } -__global__ void ftest(int diagFlag, int p, int rows, int colsx, int colsy, - int rCols, int unrCols, real * obs, int obsDim, - real * rCoeffs, int rCoeffsDim, real * unrCoeffs, int unrCoeffsDim, - real * rdata, int rdataDim, real * unrdata, int unrdataDim, - real * dfStats) // real * dpValues) +// The fast form for every (effect jj of the tile, cause i): from +// ll (lag x lag per cause), c = Q'L (b*m x p*lag), le = L'e (p*lag x b), +// qe = Q'e (m x b), l2 and rss. Writes gain, cond, left (b x p, cause +// fastest) and redo = 1 where the pair must be recomputed explicitly: its +// score cond / left * rss / (rss - gain) is not below limit (a NaN score, +// a constant series, has no test either way). Pairs of a series with +// itself (self set, cause i == effect j0 + jj) are NaN with no recompute. +__global__ void grangerFast(int b, int p, int lag, int j0, int self, + const real * ll, const real * c, const real * le, + const real * qe, const real * l2, const real * rss, + real limit, real * gain, real * cond, real * left, + int * redo) { - int - j = blockIdx.x * THREADSPERDIM + threadIdx.x, - i = blockIdx.y * THREADSPERDIM + threadIdx.y, - idx = i*colsx + j, k, m; - real - kobs, fp = (real) p, frows = (real) rows, - rSsq, unrSsq, - rEst, unrEst, - score = real(0), - * tObs, * tRCoeffs, * tUnrCoeffs, - * tRdata, * tUnrdata; - - if((i >= colsy) || (j >= colsx)) return; - if((!diagFlag) && (i == j)) { - dfStats[idx] = real(0); - // dpValues[idx] = real(0); + size_t e = blockIdx.x * (size_t) blockDim.x + threadIdx.x; + if (e >= (size_t) b * p) return; + int i = e % p, jj = e / p, m = lag + 1, ldc = b * m; + if (self && i == j0 + jj) { + gain[e] = REAL_NAN; cond[e] = REAL_NAN; left[e] = REAL_NAN; redo[e] = 0; return; } - - tObs = obs + (i*colsx+j)*obsDim; - - tRCoeffs = rCoeffs + i*rCoeffsDim; - tRdata = rdata + i*rdataDim; - - tUnrCoeffs = unrCoeffs + (i*colsx+j)*unrCoeffsDim; - tUnrdata = unrdata + (i*colsx+j)*unrdataDim; - - rSsq = unrSsq = real(0); - for(k = 0; k < rows; k++) { - unrEst = rEst = real(0); - kobs = tObs[k]; - for(m = 0; m < rCols; m++) - rEst += tRCoeffs[m] * tRdata[k+m*rows]; - for(m = 0; m < unrCols; m++) - unrEst += tUnrCoeffs[m] * tUnrdata[k+m*rows]; - rSsq += (kobs - rEst) * (kobs - rEst); - unrSsq += (kobs - unrEst) * (kobs - unrEst); - + real G[GRANGER_MAXLAG * GRANGER_MAXLAG], g[GRANGER_MAXLAG]; + const real * cb = c + (size_t) jj * m + (size_t) i * lag * ldc; // C[a, k] at cb[a + k*ldc] + const real * qeb = qe + (size_t) jj * m; + for (int k = 0; k < lag; k++) { + for (int l = 0; l < lag; l++) { + real s = ll[k + l * lag + (size_t) i * lag * lag]; + for (int a = 0; a < m; a++) s -= cb[a + (size_t) k * ldc] * cb[a + (size_t) l * ldc]; + G[k + l * lag] = s; + } + real s = le[(size_t) i * lag + k + (size_t) jj * (p * lag)]; + for (int a = 0; a < m; a++) s -= cb[a + (size_t) k * ldc] * qeb[a]; + g[k] = s; } - score = ((rSsq - unrSsq)*(frows-real(2)*fp-real(1))) / (fp*unrSsq); - - if(!isfinite(score)) - score = real(0); + real cnd, lft; + real gn = grangerSolvePair(lag, G, g, l2 + (size_t) i * lag, &cnd, &lft); + gain[e] = gn; cond[e] = cnd; left[e] = lft; + real r = rss[jj], amp = r / (r - (gn > real(0) ? gn : real(0))); + if (!(amp > real(0))) amp = REAL_INF; + real score = lft > real(0) ? cnd / lft * amp : REAL_INF; + redo[e] = (!(score < limit) && !isnan(score)) ? 1 : 0; +} - dfStats[idx] = score; +// The explicit form for nred flagged pairs (effect idxJ[r] of the tile, +// cause idxI[r]): M = L_i - Q_jj (Q_jj' L_i) formed in two passes over the +// T rows, then G = M'M and g = M'e, and the same solve. One thread per +// pair. +__global__ void grangerExplicit(int T, int p, int lag, int nred, + const int * idxJ, const int * idxI, + const real * q, const real * l, const real * e, + const real * l2, real * gain, real * cond, + real * left) +{ + int r = blockIdx.x * blockDim.x + threadIdx.x; + if (r >= nred) return; + int jj = idxJ[r], i = idxI[r], m = lag + 1; + const real * qb = q + (size_t) jj * m * T; + const real * lb = l + (size_t) i * lag * T; + const real * eb = e + (size_t) jj * T; + real cc[GRANGER_MAXM * GRANGER_MAXLAG]; // Q' L_i, m x lag + for (int k = 0; k < m * lag; k++) cc[k] = real(0); + for (int t = 0; t < T; t++) + for (int a = 0; a < m; a++) { + real qa = qb[t + (size_t) a * T]; + for (int k = 0; k < lag; k++) cc[a + k * m] += qa * lb[t + (size_t) k * T]; + } + real G[GRANGER_MAXLAG * GRANGER_MAXLAG], g[GRANGER_MAXLAG], mk[GRANGER_MAXLAG]; + for (int k = 0; k < lag * lag; k++) G[k] = real(0); + for (int k = 0; k < lag; k++) g[k] = real(0); + for (int t = 0; t < T; t++) { + for (int k = 0; k < lag; k++) { + real s = lb[t + (size_t) k * T]; + for (int a = 0; a < m; a++) s -= qb[t + (size_t) a * T] * cc[a + k * m]; + mk[k] = s; + } + real et = eb[t]; + for (int k = 0; k < lag; k++) { + g[k] += mk[k] * et; + for (int l = 0; l < lag; l++) G[k + l * lag] += mk[k] * mk[l]; + } + } + real cnd, lft; + size_t o = (size_t) jj * p + i; + gain[o] = grangerSolvePair(lag, G, g, l2 + (size_t) i * lag, &cnd, &lft); + cond[o] = cnd; left[o] = lft; } diff --git a/man/gpuGranger.Rd b/man/gpuGranger.Rd index bfe55ed..21aa862 100644 --- a/man/gpuGranger.Rd +++ b/man/gpuGranger.Rd @@ -1,70 +1,106 @@ \name{gpuGranger} \alias{gpuGranger} -\title{Perform Granger Causality Tests for Vectors on a GPU} +\title{Granger Causality Tests for Every Pair of Series on a GPU} \description{ - This function performs, with the aid of a GPU, Granger Causality Tests on - permutations of pairs of columns of the input matrices 'x' and 'y'. + Tests, for every ordered pair of time series, whether the past of one + improves the prediction of the other beyond the other's own past: the + F test of two least-squares fits that \pkg{lmtest}'s + \code{grangertest} and statsmodels' \code{grangercausalitytests} + (\code{ssr_ftest}) report. Every pair at once, on the GPU, in double + precision. } \usage{ gpuGranger(x, y=NULL, lag, - precision = getOption("gputools.precision", "single")) + precision = getOption("gputools.precision", "single")) } \arguments{ - \item{x}{a matrix of floating point values. Each column represents - a sequence of observations for a single random variable.} - \item{y}{an optional matrix of floating point values. Each column - represents a sequence of observations for a single random variable.} - \item{lag}{a positive integer by which to offset the sequence of - observations to calculate the coefficient for Granger causality.} - \item{precision}{\code{"single"} or \code{"double"}: the floating point precision - the GPU computes in. Single precision is the default, set with - \code{options(gputools.precision = "double")}; it runs many times - faster on most consumer GPUs. Results are returned as doubles either way.} + \item{x}{a numeric matrix, one series per column in time order down the + rows, with no missing values. Its columns are the causes, and with + \code{y} missing the effects too.} + \item{y}{an optional numeric matrix of effects, the same times in the + same rows.} + \item{lag}{how many past values of each series the models use: a whole + number from 1 to 16.} + \item{precision}{accepted for the common interface, but the Granger + tests always compute in double: in single precision the F statistics + of long random walks were off by more than 10 percent.} } -\value{a list of two single precision floating point matrices both of the same - dimension. The two matrices are fStatistics and pValues. The - fStatistics matrix holds the F-statistics from the Granger causality tests. - Each element of the pValues matrix is the p-value for the corresponding - element of the fStatistics matrix. +\details{ + With \eqn{T = n - lag} usable times, the restricted model of effect + \eqn{j} regresses it on an intercept and its own \eqn{lag} past values, + the unrestricted one adds cause \eqn{i}'s \eqn{lag} past values, and + \deqn{F = \frac{(RSS_r - RSS_u) / lag}{RSS_u / (T - 2\,lag - 1)}} + on \eqn{(lag, T - 2\,lag - 1)} degrees of freedom. - Both matrices have dimnames \code{cause} (rows) and \code{effect} - (columns), taken from the column names of \code{x} and \code{y}. + By the Frisch-Waugh theorem \eqn{RSS_r - RSS_u = g' G^{-1} g} where, + after the effect's restricted design is projected out of everything + (through an orthonormal basis of it: Cholesky-QR2, or Householder QR for + designs too ill-conditioned for it), \eqn{M} holds the cause's lags with + the effect's own past removed, \eqn{g = M'e} is their product with the + effect's residual and \eqn{G = M'M}. Each effect's basis is built once, + and every cause's test against it comes from a few matrix products. For + most pairs \eqn{G} and \eqn{g} are taken from the small products + \eqn{L'L - C'C} and \eqn{L'e - C'(Q'e)} with \eqn{C = Q'L}; that form + subtracts numbers of the size of the cause's lags to leave what is + beside the effect's own past, so where its rounding error could matter + (judged pair by pair from the conditioning of \eqn{G}, how much of the + cause is left, and the size of F) the pair is recomputed with \eqn{M} + formed explicitly. Each series is first centred and scaled, which + changes no F. In double precision the F statistics agree with + \code{lm()} to about 1e-8. - If y is NULL, the test is run on permutations of pairs of columns of x. To - find the Granger causality F-statistic estimating the answer to - "Does variable x[ ,i] Granger-cause variable x[ ,j]?", look at - fStatistics[i, j] and pValues[i, j]. The diagonal, a variable against - itself, is \code{NA}. + A test is \code{NA}, with a warning, when it cannot be resolved in + double precision: when the effect's restricted design (intercept and + lags, each at unit length) has condition number above 1e10, or a + cause's \eqn{G} above 1e12, or when less than 1e-20 of a cause's lags + (in squared length) is left beside the effect's own past, a cause that + is, to rounding, a linear function of it. (The condition numbers are + the estimates \eqn{\|A\|_F \|A^{-1}\|_F}, within a factor of the number + of columns of the usual ones.) White noise, near-unit-root series, + random walks with or without drift and trends stay far below these + limits; series integrated two or three times can exceed them at long + lags, and should be differenced first, as Granger tests assume + stationary series anyway. Exactly periodic series (a pure sinusoid) + have collinear lags too. - If y is not NULL, the test is run on permutations of pairs (x[ ,i], y[ ,j]). - To find the Granger causality F-statistic estimating the answer to - "Does variable x[ ,i] Granger-cause variable y[ ,j]?", look at - fStatistics[i, j] and pValues[i, j]. + A series is not tested against itself (its lags are already in the + restricted model), and a constant series has no tests: \code{NA}, with + no warning. } -\examples{ -\dontshow{if (gpuAvailable()) withAutoprint(\{ # examplesIf} -# permutations of pairs of cols of just x -numRandVars <- 5 -numSamples <- 20 -randVarSequences <- matrix(runif(numRandVars*numSamples), numSamples, - numRandVars) -gpuGranger(randVarSequences, lag = 5) +\value{a list of two numeric matrices, \code{fStatistics} and + \code{pValues}, each with a row per column of \code{x} and a column per + column of \code{y} (of \code{x} when \code{y} is missing), with dimnames + \code{cause} (rows) and \code{effect} (columns) taken from the column + names of \code{x} and \code{y}. -# pairs of cols, one from x and one from y -numXRandVars <- 5 -numXSamples <- 20 -x <- matrix(runif(numXRandVars*numXSamples), numXSamples, numXRandVars) + \code{fStatistics[i, j]} answers "does \code{x[, i]} Granger-cause + \code{y[, j]}?" (or \code{x[, j]}), and \code{pValues[i, j]} is its + p-value from the F distribution on \eqn{(lag, n - 3\,lag - 1)} degrees + of freedom. With \code{y} missing the diagonal, a series against itself, + is \code{NA}. +} -numYRandVars <- 3 -numYSamples <- 20 -y <- matrix(runif(numYRandVars*numYSamples), numYSamples, numYRandVars) +\references{ + Granger, C. W. J. (1969). Investigating causal relations by econometric + models and cross-spectral methods. \emph{Econometrica} 37(3), 424-438. +} + +\examples{ +\dontshow{if (gpuAvailable()) withAutoprint(\{ # examplesIf} +set.seed(1) +x <- matrix(rnorm(300 * 4), 300, 4) +x[2:300, 3] <- x[2:300, 3] + 0.8 * x[1:299, 1] # series 1 drives series 3 +res <- gpuGranger(x, lag = 2) +which(res$pValues < 1e-6, arr.ind = TRUE) # cause 1, effect 3 -result <- gpuGranger(x, y, lag = 5) -print(result) +# causes in x, effects in y +y <- matrix(rnorm(300 * 2), 300, 2) +y[2:300, 2] <- y[2:300, 2] + 0.5 * x[1:299, 4] +gpuGranger(x, y, lag = 1)$pValues \dontshow{\}) # examplesIf} } diff --git a/src/cudaUtils.cpp b/src/cudaUtils.cpp index a26d7df..ffe5d62 100644 --- a/src/cudaUtils.cpp +++ b/src/cudaUtils.cpp @@ -41,7 +41,9 @@ const std::vector & kernelFiles() { "euclidean_kernel_same", "maximum_kernel_same", "manhattan_kernel_same", "canberra_kernel_same", "binary_kernel_same", "minkowski_kernel_same" } }, - { "granger", { "getRestricted", "getUnrestricted", "ftest" } }, + { "granger", { "grangerDesign", "grangerFactor", "grangerApplyRinv", + "grangerRss", "grangerColumnSq", "grangerFast", + "grangerExplicit" } }, { "hcluster", { "complete_kernel", "wpgma_kernel", "average_kernel", "median_kernel", "centroid_kernel", "flexible_group_kernel", "flexible_kernel", diff --git a/src/gpublas.h b/src/gpublas.h index 89258a7..256d816 100644 --- a/src/gpublas.h +++ b/src/gpublas.h @@ -145,6 +145,19 @@ inline void trsmLeftUpper(int m, int n, float alpha, const float * a, CUBLAS_DIAG_NON_UNIT, m, n, &alpha, a, lda, b, ldb)); } +// b independent products c_k = alpha op(a_k) op(b_k) + beta c_k, the k-th +// operands strideA, strideB, strideC elements after the first +inline void gemmStridedBatched(char transa, char transb, int m, int n, int k, + float alpha, const float * a, int lda, + long long strideA, const float * b, int ldb, + long long strideB, float beta, float * c, + int ldc, long long strideC, int batch) +{ + CUBLAS_CALL(cublasSgemmStridedBatched(gpuBlasHandle(), op(transa), + op(transb), m, n, k, &alpha, a, lda, strideA, b, ldb, strideB, &beta, c, + ldc, strideC, batch)); +} + // the same in double precision inline void copy(int n, const double * x, int incx, double * y, int incy) @@ -212,6 +225,17 @@ inline void gemm(char transa, char transb, int m, int n, int k, double alpha, &alpha, a, lda, b, ldb, &beta, c, ldc)); } +inline void gemmStridedBatched(char transa, char transb, int m, int n, int k, + double alpha, const double * a, int lda, + long long strideA, const double * b, int ldb, + long long strideB, double beta, double * c, + int ldc, long long strideC, int batch) +{ + CUBLAS_CALL(cublasDgemmStridedBatched(gpuBlasHandle(), op(transa), + op(transb), m, n, k, &alpha, a, lda, strideA, b, ldb, strideB, &beta, c, + ldc, strideC, batch)); +} + inline void syrk(char uplo, char trans, int n, int k, double alpha, const double * a, int lda, double beta, double * c, int ldc) { diff --git a/src/granger.cpp b/src/granger.cpp index 179bfa3..c110cf6 100644 --- a/src/granger.cpp +++ b/src/granger.cpp @@ -1,381 +1,365 @@ -#include -#include -#include -#include +// Granger causality F tests for every ordered pair of series. +// +// For each ordered pair (cause i, effect j), with T = rows - lag usable +// times: the restricted model regresses effect j on an intercept and its +// own lag past values, the unrestricted one adds cause i's lag past +// values, and F = ((RSS_r - RSS_u) / lag) / (RSS_u / (T - 2 lag - 1)). +// +// By the Frisch-Waugh theorem RSS_r - RSS_u = g' G^{-1} g where, after +// projecting the effect's restricted design out of everything through an +// orthonormal basis Q of it, M = L_i - Q Q' L_i holds the cause's lags with +// the effect's own past removed, G = M'M and g = M'e for the effect's +// residual e. Each effect's Q is built once (Cholesky-QR2 on the GPU, +// Householder on the host for designs too ill-conditioned for it), then +// every cause's test against it comes from a few matrix products: the fast +// form G = L'L - C'C, g = L'e - C'(Q'e) with C = Q'L, which the kernel +// checks pair by pair and falls back from to M formed explicitly where the +// cancellation cannot be trusted (see FAST_SCORE_LIMIT). Everything is +// computed in double; each series is centred and scaled first, which +// changes no F. +// +// The method, its limits and the measurements behind them are those of +// the Python package widestats (granger.py), which this follows. + +#include +#include +#include +#include #include "Rmath.h" -#include "cuda_runtime_api.h" +#include "cuda_runtime.h" -#include "cuseful.h" #include "cudaUtils.h" +#include "gpublas.h" +#include "cuseful.h" #include "granger.h" -#define max(a, b) ((a > b)?a:b) - -#define THREADSPERDIM 16 - -#define FALSE 0 -#define TRUE !FALSE +namespace { + +const double NaN = std::numeric_limits::quiet_NaN(); +const double Inf = std::numeric_limits::infinity(); + +// Condition numbers (Frobenius estimates ||A|| ||A^{-1}||, at least the +// 2-norm ones and at most the number of columns times them) above which a +// test is NaN: of the effect's restricted design, and of a cause's G; and +// the fraction of a cause's lags (squared length) left beside the effect's +// own past below which nothing is left to test. Measured against least +// squares over 1,152 pairs (white noise, AR(0.9999), random walks with and +// without drift, trends, I(2) with and without drift, I(3), mixed; n +// 300-30,000; lag 1-8): realistic series stay below 5e4 and 3e5 with +// errors under 1e-8; within these limits every pair was within 1e-4; +// beyond them errors reach 100%. +const double EFFECT_LIMIT = 1e10, CAUSE_LIMIT = 1e12, LEFT_LIMIT = 1e-20; +// Cholesky-QR2 for designs better conditioned than this (as accurate as +// Householder there, measured); Householder QR for the rest. +const double CHOLESKY_QR_LIMIT = 1e7; +// The fast form's error relative to the explicit one grows with +// cond(G) / left and with RSS_r / RSS_u; over 1,501 resolved pairs it was +// at most about 20 eps times that product, so pairs with the product at or +// above this are recomputed from M (error below 2.2e-10 measured). +const double FAST_SCORE_LIMIT = 1e6; +const int MAXLAG = 16; // GRANGER_MAXLAG in granger.cu template -void getPValues(int rows, int cols, const T * fStats, int p, int embedRows, - T * pValues) +T * upload(const std::vector & h) { - T fscore = T(0); - - for(int i = 0; i < rows; i++) { - for(int j = 0; j < cols; j++) { - fscore = fStats[i + j * rows]; - // the upper tail directly: 1 - pf(lower tail) rounds small p to 0 - pValues[i + j * rows] = (T) pf((double) fscore, (double) p, - (double) embedRows - 2.0 * p - 1.0, - 0, 0); - } - } + T * d = gpuAlloc(h.size()); + CUDART_CALL(cudaMemcpy(d, h.data(), h.size() * sizeof(T), + cudaMemcpyHostToDevice)); + return d; } template -void granger(int rows, int cols, const T * y, int p, - T * fStats, T * pValues) +void download(const T * d, std::vector & h) { - if(cols < 2) { - fatal("The Granger test needs at least 2 variables.\n"); - return; - } - int - i, j, k, t = p+1, - fbytes = sizeof(T), - embedRows = rows-p, embedCols = t*2; - T - * Y, * rQ, * rR, - * unrQ, * unrR, - * restricted, * unrestricted, - * rdata, * unrdata, - * dfStats; // * dpValues; - size_t - size = cols*cols*fbytes, partSize = embedRows*size; - - cudaMalloc((void **)&Y, embedCols*partSize); - - cudaMalloc((void **)&rQ, t*embedRows*cols*fbytes); - cudaMalloc((void **)&rR, t*t*cols*fbytes); - cudaMalloc((void **)&rdata, t*embedRows*cols*fbytes); - cudaMalloc((void **)&unrdata, (embedCols-1)*partSize); - cudaMalloc((void **)&restricted, t*cols*fbytes); - if( hasCudaError("granger: line 267: gpu memory allocation") ) return; + CUDART_CALL(cudaMemcpy(h.data(), d, h.size() * sizeof(T), + cudaMemcpyDeviceToHost)); +} - int - Ydim = embedCols * embedRows, - rQdim = t * embedRows, rRdim = t * t, - rdataDim = t*embedRows, restrictedDim = t, - unrQdim = (embedCols-1) * embedRows, - unrRdim = (embedCols-1) * (embedCols-1), - unrestrictedDim = embedCols-1, unrdataDim = (embedCols-1)*embedRows; - T - * ypos, * rdataPos, * unrdataPos, - * evenCols; - int - skip = 2*embedRows, colBytes = embedRows*fbytes; - const T - * vectA, * vectB; +dim3 gridFor(size_t n, int threads) +{ + return dim3((unsigned) ((n + threads - 1) / threads)); +} - for(i = 0; i < cols; i++) { - rdataPos = rdata+i*rdataDim; - evenCols = rdataPos+embedRows; - vectA = y+i*rows; - for(j = 0; j < cols; j++) { - if(i == j) continue; - - ypos = Y+(i*cols+j)*Ydim; - unrdataPos = unrdata+(i*cols+j)*unrdataDim; +// Centres and scales each column of a (rows x cols) to unit length; a +// constant column becomes NaN. +void scale(int rows, int cols, const double * a, std::vector & out) +{ + out.assign((size_t) rows * cols, 0.0); + for (int j = 0; j < cols; j++) { + const double * col = a + (size_t) j * rows; + double mean = 0; + for (int t = 0; t < rows; t++) mean += col[t]; + mean /= rows; + double ss = 0; + for (int t = 0; t < rows; t++) ss += (col[t] - mean) * (col[t] - mean); + double s = std::sqrt(ss); + for (int t = 0; t < rows; t++) + out[(size_t) j * rows + t] = (col[t] - mean) / s; + } +} - vectB = y+j*rows; +// The lag matrix of the scaled series s (rows x cols): T x cols*lag, +// column j*lag + k the value k + 1 steps back of series j. +void lagged(int rows, int cols, int lag, const std::vector & s, + std::vector & out) +{ + int T = rows - lag; + out.assign((size_t) T * cols * lag, 0.0); + for (int j = 0; j < cols; j++) + for (int k = 0; k < lag; k++) { + const double * src = s.data() + (size_t) j * rows + lag - k - 1; + double * dst = out.data() + ((size_t) j * lag + k) * T; + for (int t = 0; t < T; t++) dst[t] = src[t]; + } +} - for(k = 0; k < p+1; k++) { // produce t subcols - cudaMemcpy(ypos+k*skip, vectA+(p-k), embedRows*fbytes, - cudaMemcpyHostToDevice); - cudaMemcpy(ypos+k*skip+embedRows, vectB+(p-k), - embedRows*fbytes, cudaMemcpyHostToDevice); - } - cudaMemcpy(unrdataPos+embedRows, ypos+skip, - (embedCols-2)*embedRows*fbytes, cudaMemcpyDeviceToDevice); +// Householder QR of the T x m column-major matrix d (on the host), for the +// designs Cholesky-QR2 cannot handle: Q (T x m) into q, and the Frobenius +// condition estimate of R returned (infinite if R is singular or not +// finite). +double householder(int T, int m, std::vector & d, + std::vector & q) +{ + std::vector v((size_t) T * m), tau(m); + for (int j = 0; j < m; j++) { + double * col = d.data() + (size_t) j * T; + double norm = 0; + for (int t = j; t < T; t++) norm += col[t] * col[t]; + norm = std::sqrt(norm); + double alpha = col[j] >= 0 ? -norm : norm; + double * vj = v.data() + (size_t) j * T; + for (int t = 0; t < T; t++) vj[t] = t < j ? 0 : col[t]; + vj[j] -= alpha; + double vv = 0; + for (int t = j; t < T; t++) vv += vj[t] * vj[t]; + tau[j] = vv > 0 ? 2 / vv : 0; + // apply H_j = I - tau v v' to the remaining columns + for (int c = j; c < m; c++) { + double * cc = d.data() + (size_t) c * T, s = 0; + for (int t = j; t < T; t++) s += vj[t] * cc[t]; + s *= tau[j]; + for (int t = j; t < T; t++) cc[t] -= s * vj[t]; } - // build restricted data from last set of unrestricted data - // only need one per column, not one for each pairing - for(k = 0; k < embedCols-2; k+=2) { - cudaMemcpy(evenCols+(k*embedRows)/2, unrdataPos+(1+k)*embedRows, - colBytes, cudaMemcpyDeviceToDevice); + } + // R is now the upper triangle of d; Q = H_0 ... H_{m-1} [I; 0] + q.assign((size_t) T * m, 0.0); + for (int c = 0; c < m; c++) q[(size_t) c * T + c] = 1; + for (int j = m - 1; j >= 0; j--) { + const double * vj = v.data() + (size_t) j * T; + for (int c = 0; c < m; c++) { + double * qc = q.data() + (size_t) c * T, s = 0; + for (int t = j; t < T; t++) s += vj[t] * qc[t]; + s *= tau[j]; + for (int t = j; t < T; t++) qc[t] -= s * vj[t]; } } - if( hasCudaError("granger : mem copy from host to device") ) return; - - int numBlocks = cols / THREADSPERDIM; - if(numBlocks * THREADSPERDIM < cols) numBlocks++; - - dim3 - dimRGrid(numBlocks), - dimRBlock(THREADSPERDIM), - dimUnrGrid(numBlocks, numBlocks), - dimUnrBlock(THREADSPERDIM, THREADSPERDIM); - - void * restArgs[] = { - &cols, &cols, - &embedRows, - &t, - &rdata, - &rdataDim, - &Y, - &Ydim, - &rQ, - &rQdim, - &rR, - &rRdim, - &restricted, - &restrictedDim - }; - cudaLaunch("getRestricted", restArgs, dimRGrid, dimRBlock, precisionOf()); - if( hasCudaError("granger: getRestricted kernel execution") ) return; - - cudaFree(rQ); - cudaFree(rR); - - cudaMalloc((void **)&unrQ, (embedCols-1)*partSize); - cudaMalloc((void **)&unrR, (embedCols-1)*(embedCols-1)*size); - cudaMalloc((void **)&unrestricted, (embedCols-1)*size); - if( hasCudaError("granger: line 336: attemped gpu memory allocation") ) - return; - - int unrestT = embedCols - 1; // the kernels take an int - void * unrestArgs[] = { - &cols, &cols, - &embedRows, - &unrestT, - &unrdata, - &unrdataDim, - &Y, - &Ydim, - &unrQ, - &unrQdim, - &unrR, - &unrRdim, - &unrestricted, - &unrestrictedDim - }; - cudaLaunch("getUnrestricted", unrestArgs, dimUnrGrid, dimUnrBlock, precisionOf()); - if( hasCudaError("granger : getUnRestricted kernel execution") ) return; - - cudaFree(unrQ); - cudaFree(unrR); - - size_t resultSize = cols*cols*fbytes; - cudaMalloc((void **)&dfStats, resultSize); - // cudaMalloc((void **)&dpValues, resultSize); - if( hasCudaError("granger: line 350: gpu memory allocation") ) return; - - int diagFlag = FALSE; - void * ftestArgs[] = { - &diagFlag, - &p, - &embedRows, - &cols, &cols, - &t, &unrestT, - &Y, &Ydim, - &restricted, &restrictedDim, - &unrestricted, &unrestrictedDim, - &rdata, &rdataDim, - &unrdata, &unrdataDim, - &dfStats - }; - cudaLaunch("ftest", ftestArgs, dimUnrGrid, dimUnrBlock, precisionOf()); - if( hasCudaError("granger : ftest kernel execution") ) return; - - cudaMemcpy(fStats, dfStats, resultSize, cudaMemcpyDeviceToHost); - // cudaMemcpy(pValues, dpValues, resultSize, cudaMemcpyDeviceToHost); - if( hasCudaError("granger : mem copy device to host") ) return; - - getPValues(cols, cols, fStats, p, embedRows, pValues); - - cudaFree(Y); - cudaFree(restricted); - cudaFree(unrestricted); - cudaFree(rdata); - cudaFree(unrdata); - cudaFree(dfStats); - // cudaFree(dpValues); + // ||R||_F ||R^{-1}||_F by back substitution against the identity + double rf = 0, rif = 0; + for (int j = 0; j < m; j++) + for (int i = 0; i <= j; i++) rf += d[i + (size_t) j * T] * d[i + (size_t) j * T]; + for (int c = 0; c < m; c++) { + std::vector x(m, 0.0); + for (int i = m - 1; i >= 0; i--) { + double s = i == c ? 1 : 0; + for (int k = i + 1; k < m; k++) s -= d[i + (size_t) k * T] * x[k]; + x[i] = s / d[i + (size_t) i * T]; + rif += x[i] * x[i]; + } + } + double cond = std::sqrt(rf) * std::sqrt(rif); + return std::isfinite(cond) ? cond : Inf; } -template -void grangerxy(int rows, int colsx, const T * x, int colsy, - const T * y, int p, T * fStats, T * pValues) -{ +} // namespace - if((p < 0) || (rows < 1) || (colsx < 1) || (colsy < 1)) { - fatal("The Granger XY test needs at least a pair variables.\n"); - return; +void granger(int rows, int colsx, const double * x, int colsy, + const double * y, int lag, double * fStats, double * pValues, + int * unresolved) +{ + const bool self = y == NULL || y == x; + if (self) { y = x; colsy = colsx; } + if (lag < 1 || lag > MAXLAG) + error("gputools: lag must be between 1 and %d", MAXLAG); + const int T = rows - lag, m = lag + 1, df2 = T - 2 * lag - 1, p = colsx, + q = colsy, pL = p * lag; + if (df2 < 1) + error("gputools: the series are too short for lag %d: need at least %d " + "times", lag, 3 * lag + 2); + const Precision prec = Precision::Double; + const int threads = 256; + + // the scaled series' lags and the effects' targets, once, on the device + std::vector xs, ys, hL, hLY, hTarget((size_t) T * q); + scale(rows, p, x, xs); + lagged(rows, p, lag, xs, hL); + if (self) { hLY = hL; } else { scale(rows, q, y, ys); lagged(rows, q, lag, ys, hLY); } + const std::vector & yscaled = self ? xs : ys; + for (int j = 0; j < q; j++) + for (int t = 0; t < T; t++) + hTarget[(size_t) j * T + t] = yscaled[(size_t) j * rows + lag + t]; + double * dL = upload(hL), * dLY = self ? dL : upload(hLY), + * dTarget = upload(hTarget); + // each cause's lag columns' squared lengths and cross-products L_i'L_i + double * dL2 = gpuAlloc(pL), * dLL = gpuAlloc((size_t) p * lag * lag); + { + int ncol = pL, TT = T; + void * args[] = { &TT, &ncol, &dL, &dL2 }; + cudaLaunch("grangerColumnSq", args, gridFor(pL, threads), dim3(threads), prec); + blas::gemmStridedBatched('T', 'N', lag, lag, T, 1.0, dL, T, (long long) lag * T, + dL, T, (long long) lag * T, 0.0, dLL, lag, + (long long) lag * lag, p); } - int - i, j, k, t = p+1, - fbytes = sizeof(T), - embedRows = rows-p, embedCols = t*2; - T - * Y, * rQ, * rR, - * unrQ, * unrR, - * restricted, * unrestricted, - * rdata, * unrdata, - * dfStats; // * dpValues; - size_t - size = colsx*colsy*fbytes, partSize = embedRows*size; - - cudaMalloc((void **)&Y, embedCols*partSize); - - cudaMalloc((void **)&rQ, t*embedRows*colsy*fbytes); - cudaMalloc((void **)&rR, t*t*colsy*fbytes); - cudaMalloc((void **)&rdata, t*embedRows*colsy*fbytes); - cudaMalloc((void **)&restricted, t*colsy*fbytes); - - cudaMalloc((void **)&unrQ, (embedCols-1)*partSize); - cudaMalloc((void **)&unrR, (embedCols-1)*(embedCols-1)*size); - cudaMalloc((void **)&unrestricted, (embedCols-1)*size); - cudaMalloc((void **)&unrdata, (embedCols-1)*partSize); - checkCudaError("grangerxy : attemped gpu memory allocation"); - int - Ydim = embedCols * embedRows, - rQdim = t * embedRows, rRdim = t * t, - rdataDim = t*embedRows, restrictedDim = t, - unrQdim = (embedCols-1) * embedRows, - unrRdim = (embedCols-1) * (embedCols-1), - unrestrictedDim = embedCols-1, unrdataDim = (embedCols-1)*embedRows; - T - * ypos, * rdataPos, * unrdataPos; - - int - skip = 2*embedRows, colBytes = embedRows*fbytes; - const T * vectA, * vectB; - T * evenCols; - - for(i = 0; i < colsy; i++) { - rdataPos = rdata+i*rdataDim; - evenCols = rdataPos+embedRows; - vectA = y+i*rows; - for(j = 0; j < colsx; j++) { - ypos = Y+(i*colsx+j)*Ydim; - unrdataPos = unrdata+(i*colsx+j)*unrdataDim; - - vectB = x+j*rows; - - for(k = 0; k < p+1; k++) { // produce t subcols - cudaMemcpy(ypos+k*skip, vectA+(p-k), embedRows*fbytes, - cudaMemcpyHostToDevice); - cudaMemcpy(ypos+k*skip+embedRows, vectB+(p-k), - embedRows*fbytes, cudaMemcpyHostToDevice); + // effects a tile at a time: the fast form holds about six T x m per + // effect (design, Q, intermediates) and C and G + size_t perEffect = (size_t) 6 * T * m + (size_t) pL * (m + lag); + int tile = (int) std::max(1, std::min(q, ((size_t) 1 << 24) / perEffect)); + double * dD = gpuAlloc((size_t) T * tile * m), + * dQ = gpuAlloc((size_t) T * tile * m), + * dGram = gpuAlloc((size_t) tile * m * m), + * dRinv = gpuAlloc((size_t) tile * m * m), + * dEffectCond = gpuAlloc(tile), + * dE = gpuAlloc((size_t) T * tile), + * dQy = gpuAlloc((size_t) tile * m), + * dQe = gpuAlloc((size_t) tile * m), + * dRss = gpuAlloc(tile), + * dC = gpuAlloc((size_t) tile * m * pL), + * dLe = gpuAlloc((size_t) pL * tile), + * dGain = gpuAlloc((size_t) tile * p), + * dCond = gpuAlloc((size_t) tile * p), + * dLeft = gpuAlloc((size_t) tile * p); + int * dBad = gpuAlloc(tile), * dRedo = gpuAlloc((size_t) tile * p), + * dIdxJ = gpuAlloc((size_t) tile * p), * dIdxI = gpuAlloc((size_t) tile * p); + + std::vector gain((size_t) p * q), causeCond((size_t) p * q), + left((size_t) p * q), rssR(q), effectCond(q); + std::vector hGain((size_t) tile * p), hCond((size_t) tile * p), + hLeft((size_t) tile * p), hEffCond(tile), hRss(tile); + std::vector hBad(tile), hRedo((size_t) tile * p), hIdxJ, hIdxI; + + for (int j0 = 0; j0 < q; j0 += tile) { + int b = std::min(tile, q - j0); + size_t bm = (size_t) b * m; + // the restricted designs, orthonormalised by Cholesky-QR2 + { + int TT = T, lg = lag, jj0 = j0, bb = b; + void * args[] = { &TT, &lg, &jj0, &bb, &dLY, &dD }; + cudaLaunch("grangerDesign", args, gridFor((size_t) T * bm, threads), dim3(threads), prec); + } + CUDART_CALL(cudaMemcpy(dQ, dD, (size_t) T * bm * sizeof(double), cudaMemcpyDeviceToDevice)); + CUDART_CALL(cudaMemset(dBad, 0, b * sizeof(int))); + for (int pass = 0; pass < 2; pass++) { + blas::gemmStridedBatched('T', 'N', m, m, T, 1.0, dQ, T, (long long) m * T, + dQ, T, (long long) m * T, 0.0, dGram, m, + (long long) m * m, b); + int bb = b, mm = m, writeCond = pass == 0; + double limit = CHOLESKY_QR_LIMIT; + void * fargs[] = { &bb, &mm, &dGram, &limit, &writeCond, &dRinv, &dEffectCond, &dBad }; + cudaLaunch("grangerFactor", fargs, gridFor(b, threads), dim3(threads), prec); + int TT = T; + void * aargs[] = { &TT, &mm, &bb, &dRinv, &dQ }; + cudaLaunch("grangerApplyRinv", aargs, gridFor((size_t) T * b, threads), dim3(threads), prec); + } + // Householder on the host for the designs Cholesky-QR2 could not take + download(dBad, hBad); + hBad.resize(b); download(dEffectCond, hEffCond); hEffCond.resize(b); + for (int jj = 0; jj < b; jj++) { + if (!hBad[jj]) continue; + std::vector d((size_t) T * m), qh; + CUDART_CALL(cudaMemcpy(d.data(), dD + (size_t) jj * m * T, d.size() * sizeof(double), + cudaMemcpyDeviceToHost)); + hEffCond[jj] = householder(T, m, d, qh); + CUDART_CALL(cudaMemcpy(dQ + (size_t) jj * m * T, qh.data(), qh.size() * sizeof(double), + cudaMemcpyHostToDevice)); + } + hBad.assign(tile, 0); hEffCond.resize(tile); + std::copy(hEffCond.begin(), hEffCond.begin() + b, effectCond.begin() + j0); + // residuals e = y - Q (Q'y), their sums of squares + CUDART_CALL(cudaMemcpy(dE, dTarget + (size_t) j0 * T, (size_t) T * b * sizeof(double), + cudaMemcpyDeviceToDevice)); + blas::gemmStridedBatched('T', 'N', m, 1, T, 1.0, dQ, T, (long long) m * T, + dE, T, (long long) T, 0.0, dQy, m, (long long) m, b); + blas::gemmStridedBatched('N', 'N', T, 1, m, -1.0, dQ, T, (long long) m * T, + dQy, m, (long long) m, 1.0, dE, T, (long long) T, b); + { + int TT = T, bb = b; + void * args[] = { &TT, &bb, &dE, &dRss }; + cudaLaunch("grangerRss", args, gridFor(b, threads), dim3(threads), prec); + } + // the fast form's ingredients: C = Q'L, L'e, Q'e + blas::gemm('T', 'N', (int) bm, pL, T, 1.0, dQ, T, dL, T, 0.0, dC, (int) bm); + blas::gemm('T', 'N', pL, b, T, 1.0, dL, T, dE, T, 0.0, dLe, pL); + blas::gemmStridedBatched('T', 'N', m, 1, T, 1.0, dQ, T, (long long) m * T, + dE, T, (long long) T, 0.0, dQe, m, (long long) m, b); + { + int bb = b, pp = p, lg = lag, jj0 = j0, sf = self; + double limit = FAST_SCORE_LIMIT; + void * args[] = { &bb, &pp, &lg, &jj0, &sf, &dLL, &dC, &dLe, &dQe, &dL2, &dRss, + &limit, &dGain, &dCond, &dLeft, &dRedo }; + cudaLaunch("grangerFast", args, gridFor((size_t) b * p, threads), dim3(threads), prec); + } + // the pairs the fast form cannot be trusted on, explicitly + hRedo.resize((size_t) b * p); download(dRedo, hRedo); + hIdxJ.clear(); hIdxI.clear(); + for (size_t e = 0; e < (size_t) b * p; e++) + if (hRedo[e]) { hIdxJ.push_back((int) (e / p)); hIdxI.push_back((int) (e % p)); } + hRedo.resize((size_t) tile * p); + int nred = (int) hIdxJ.size(); + if (nred) { + CUDART_CALL(cudaMemcpy(dIdxJ, hIdxJ.data(), nred * sizeof(int), cudaMemcpyHostToDevice)); + CUDART_CALL(cudaMemcpy(dIdxI, hIdxI.data(), nred * sizeof(int), cudaMemcpyHostToDevice)); + int TT = T, pp = p, lg = lag; + void * args[] = { &TT, &pp, &lg, &nred, &dIdxJ, &dIdxI, &dQ, &dL, &dE, &dL2, + &dGain, &dCond, &dLeft }; + cudaLaunch("grangerExplicit", args, gridFor(nred, threads), dim3(threads), prec); + } + hGain.resize((size_t) b * p); download(dGain, hGain); + hCond.resize((size_t) b * p); download(dCond, hCond); + hLeft.resize((size_t) b * p); download(dLeft, hLeft); + hRss.resize(b); download(dRss, hRss); + for (int jj = 0; jj < b; jj++) { + rssR[j0 + jj] = hRss[jj]; + for (int i = 0; i < p; i++) { + size_t src = (size_t) jj * p + i, dst = (size_t) (j0 + jj) * p + i; + gain[dst] = hGain[src]; causeCond[dst] = hCond[src]; left[dst] = hLeft[src]; } - cudaMemcpy(unrdataPos+embedRows, ypos+skip, - (embedCols-2)*embedRows*fbytes, cudaMemcpyDeviceToDevice); } - // build restricted data from last set of unrestricted data - // only need one per column, not one for each pairing - for(k = 0; k < embedCols-2; k+=2) { - cudaMemcpy(evenCols+(k*embedRows)/2, unrdataPos+(1+k)*embedRows, - colBytes, cudaMemcpyDeviceToDevice); + hGain.resize((size_t) tile * p); hCond.resize((size_t) tile * p); + hLeft.resize((size_t) tile * p); hRss.resize(tile); + } + // free before the error check: it raises (longjmps) if a kernel faulted + cudaFree(dD); cudaFree(dQ); cudaFree(dGram); cudaFree(dRinv); cudaFree(dEffectCond); + cudaFree(dE); cudaFree(dQy); cudaFree(dQe); cudaFree(dRss); cudaFree(dC); cudaFree(dLe); + cudaFree(dGain); cudaFree(dCond); cudaFree(dLeft); cudaFree(dBad); cudaFree(dRedo); + cudaFree(dIdxJ); cudaFree(dIdxI); cudaFree(dL2); cudaFree(dLL); + if (!self) cudaFree(dLY); + cudaFree(dL); cudaFree(dTarget); + if (hasCudaError("granger")) { /* the error is raised by the check */ } + + // the tests, NaN where unresolvable: a series on itself and a constant + // series have no test (NaN, not counted); the rest are counted for the + // caller's warning + int count = 0; + for (int j = 0; j < q; j++) { + bool effectWeak = std::isfinite(effectCond[j]) && effectCond[j] > EFFECT_LIMIT; + for (int i = 0; i < p; i++) { + size_t e = (size_t) j * p + i; + double gn = gain[e]; + bool weak = (effectWeak || causeCond[e] > CAUSE_LIMIT || left[e] < LEFT_LIMIT) + && (std::isfinite(gn) || std::isinf(causeCond[e])); + if (self && i == j) weak = false; + if (weak) { gn = NaN; count++; } + if (std::isnan(gn)) { fStats[e] = NaN; pValues[e] = NaN; continue; } + // rounding can take a gain of 0 just below 0, or one equal to the + // restricted RSS just above it + gn = std::min(std::max(gn, 0.0), rssR[j]); + double rssU = rssR[j] - gn; + double f = (gn / lag) / (rssU / df2); + fStats[e] = f; + // the upper tail directly: 1 - pf(lower tail) rounds small p to 0 + pValues[e] = pf(f, (double) lag, (double) df2, 0, 0); } - if( hasCudaError("grangerxy : building lagged data") ) return; } - checkCudaError("grangerxy : mem copy from host to device"); - - int - numBlocksX = colsx / THREADSPERDIM, - numBlocksY = colsy / THREADSPERDIM; - - if(numBlocksX * THREADSPERDIM < colsx) numBlocksX++; - if(numBlocksY * THREADSPERDIM < colsy) numBlocksY++; - - dim3 - dimRGrid(numBlocksY), - dimRBlock(THREADSPERDIM), - dimUnrGrid(numBlocksX, numBlocksY), - dimUnrBlock(THREADSPERDIM, THREADSPERDIM); - - void * restArgs[] = { - &colsx, &colsy, - &embedRows, - &t, - &rdata, &rdataDim, - &Y, &Ydim, - &rQ, &rQdim, - &rR, &rRdim, - &restricted, &restrictedDim - }; - cudaLaunch("getRestricted", restArgs, - dimRGrid, dimRBlock, precisionOf()); - - int unrestT = embedCols - 1; // the kernels take an int - void * unrestArgs[] = { - &colsx, &colsy, - &embedRows, - &unrestT, - &unrdata, &unrdataDim, - &Y, &Ydim, - &unrQ, &unrQdim, - &unrR, &unrRdim, - &unrestricted, &unrestrictedDim - }; - cudaLaunch("getUnrestricted", unrestArgs, - dimUnrGrid, dimUnrBlock, precisionOf()); - - checkCudaError("grangerxy : kernel execution get(Un)Restricted"); - - cudaFree(rQ); - cudaFree(unrQ); - cudaFree(rR); - cudaFree(unrR); - - size_t resultSize = colsx*colsy*fbytes; - cudaMalloc((void **)&dfStats, resultSize); - // cudaMalloc((void **)&dpValues, resultSize); - checkCudaError("grangerxy : attemped gpu memory allocation"); - - int diagFlag = TRUE; - void * ftestArgs[] = { - &diagFlag, - &p, - &embedRows, - &colsx, &colsy, - &t, &unrestT, - &Y, &Ydim, - &restricted, &restrictedDim, - &unrestricted, &unrestrictedDim, - &rdata, &rdataDim, - &unrdata, &unrdataDim, - &dfStats - }; - cudaLaunch("ftest", ftestArgs, dimUnrGrid, dimUnrBlock, precisionOf()); - checkCudaError("grangerxy : kernel execution ftest"); - - cudaMemcpy(fStats, dfStats, resultSize, cudaMemcpyDeviceToHost); - // cudaMemcpy(pValues, dpValues, resultSize, cudaMemcpyDeviceToHost); - checkCudaError("grangerxy : mem copy from device to host"); - - getPValues(colsx, colsy, fStats, p, embedRows, pValues); - - cudaFree(Y); - cudaFree(restricted); - cudaFree(unrestricted); - cudaFree(rdata); - cudaFree(unrdata); - cudaFree(dfStats); - // cudaFree(dpValues); + *unresolved = count; } - -template void getPValues(int rows, int cols, const float * fStats, int p, int embedRows, - float * pValues); -template void getPValues(int rows, int cols, const double * fStats, int p, int embedRows, - double * pValues); -template void granger(int rows, int cols, const float * y, int p, - float * fStats, float * pValues); -template void granger(int rows, int cols, const double * y, int p, - double * fStats, double * pValues); -template void grangerxy(int rows, int colsx, const float * x, int colsy, - const float * y, int p, float * fStats, float * pValues); -template void grangerxy(int rows, int colsx, const double * x, int colsy, - const double * y, int p, double * fStats, double * pValues); diff --git a/src/granger.h b/src/granger.h index fc03325..0dad62f 100644 --- a/src/granger.h +++ b/src/granger.h @@ -1,16 +1,14 @@ #ifndef _GRANGER_H_ #define _GRANGER_H_ -template -void granger(int rows, int cols, - const T * y, int p, - T * fStats, T * pValues); -template -void grangerxy(int rows, - int colsx, const T * x, - int colsy, const T * y, - int p, - T * fStats, - T * pValues); +// Granger causality F tests of every cause (column of x) on every effect +// (column of y; x itself when y is NULL), lag past values each, in double. +// fStats and pValues are colsx x colsy, column-major, [i + j*colsx] the +// test of cause i on effect j; NaN where there is no test (a series on +// itself, a constant series) or it cannot be resolved in float64, and +// unresolved counts the latter. +void granger(int rows, int colsx, const double * x, int colsy, + const double * y, int lag, double * fStats, double * pValues, + int * unresolved); #endif /* _GRANGER_H_ */ diff --git a/src/rinterface.cpp b/src/rinterface.cpp index c3e7f32..c0ce799 100644 --- a/src/rinterface.cpp +++ b/src/rinterface.cpp @@ -519,35 +519,26 @@ SEXP rankCorT(const char * kernel, SEXP x, int na, SEXP y, int nb, int size) return out; } -template -SEXP fAndP(const std::vector & f, const std::vector & p) -{ - const char * names[] = { "fStats", "pValues", "" }; +// Granger tests always compute in double (see granger.cpp) +SEXP grangerD(SEXP x, int rows, int nx, SEXP y, int ny, int lag) +{ + std::vector a = narrow(x, "x"), b; + const double * py = NULL; + if (!isNull(y)) { b = narrow(y, "y"); py = b.data(); } + else ny = nx; + std::vector f((size_t) nx * ny), p((size_t) nx * ny); + int unresolved = 0; + granger(rows, nx, a.data(), ny, py, lag, f.data(), p.data(), &unresolved); + const char * names[] = { "fStats", "pValues", "unresolved", "" }; SEXP values[] = { - PROTECT(doubles(f.data(), f.size())), PROTECT(doubles(p.data(), p.size())) + PROTECT(doubles(f.data(), f.size())), PROTECT(doubles(p.data(), p.size())), + PROTECT(ScalarInteger(unresolved)) }; SEXP out = namedList(names, values); - UNPROTECT(2); + UNPROTECT(3); return out; } -template -SEXP grangerT(SEXP x, int rows, int n, int lag) -{ - std::vector a = narrow(x, "x"), f(n * n), p(n * n); - granger(rows, n, a.data(), lag, f.data(), p.data()); - return fAndP(f, p); -} - -template -SEXP grangerXYT(SEXP x, int rows, int nx, SEXP y, int ny, int lag) -{ - std::vector - a = narrow(x, "x"), b = narrow(y, "y"), f(nx * ny), p(nx * ny); - grangerxy(rows, nx, a.data(), ny, b.data(), lag, f.data(), p.data()); - return fAndP(f, p); -} - template SEXP mutualInfoT(int bins, int order, int nsamples, SEXP x, int na, SEXP y, int nb) @@ -854,19 +845,16 @@ SEXP gpuSpearman(SEXP x, SEXP nx, SEXP y, SEXP ny, SEXP size, SEXP gpuGranger(SEXP x, SEXP rows, SEXP cols, SEXP lag, SEXP precision) { requireDevice(); - int r = asInteger(rows), n = asInteger(cols), p = asInteger(lag); - return isDouble(precision) ? grangerT(x, r, n, p) - : grangerT(x, r, n, p); + return grangerD(x, asInteger(rows), asInteger(cols), R_NilValue, 0, + asInteger(lag)); } SEXP gpuGrangerXY(SEXP x, SEXP rows, SEXP colsx, SEXP y, SEXP colsy, SEXP lag, SEXP precision) { requireDevice(); - int r = asInteger(rows), nx = asInteger(colsx), ny = asInteger(colsy), - p = asInteger(lag); - return isDouble(precision) ? grangerXYT(x, r, nx, y, ny, p) - : grangerXYT(x, r, nx, y, ny, p); + return grangerD(x, asInteger(rows), asInteger(colsx), y, asInteger(colsy), + asInteger(lag)); } SEXP gpuDist(SEXP points, SEXP numPoints, SEXP dim, SEXP method, SEXP p, diff --git a/tests/testthat/helper-gpu.R b/tests/testthat/helper-gpu.R index 15141a7..35798bf 100644 --- a/tests/testthat/helper-gpu.R +++ b/tests/testthat/helper-gpu.R @@ -9,13 +9,15 @@ skip_if_no_gpu <- function() { } # Granger causality F-test, "does cause[] Granger-cause effect[]", by lm(). -grangerRef <- function(effect, cause, lag) { +# tol is lm()'s rank tolerance: its default 1e-7 drops a cause within 1e-6 +# of the effect as aliased (F = 0), where the fits are still well defined +grangerRef <- function(effect, cause, lag, tol = 1e-7) { n <- length(effect) rows <- (lag + 1):n lags <- function(v) sapply(1:lag, function(k) v[rows - k]) y <- effect[rows] - restricted <- lm(y ~ lags(effect)) - unrestricted <- lm(y ~ lags(effect) + lags(cause)) + restricted <- lm(y ~ lags(effect), tol = tol) + unrestricted <- lm(y ~ lags(effect) + lags(cause), tol = tol) rssR <- sum(residuals(restricted)^2) rssU <- sum(residuals(unrestricted)^2) df2 <- length(rows) - 2 * lag - 1 diff --git a/tests/testthat/test-granger-mi.R b/tests/testthat/test-granger-mi.R index c12bdf8..099f213 100644 --- a/tests/testthat/test-granger-mi.R +++ b/tests/testthat/test-granger-mi.R @@ -6,42 +6,141 @@ noise <- rnorm(len) x <- cbind(effect, cause, noise) # fStatistics[i, j] answers "does column i Granger-cause column j?" +# The tests compute in double whatever the precision argument says for (pr in precisions) { - test_that(paste("gpuGranger matches an lm() F-test in", pr), { + test_that(paste("gpuGranger matches an lm() F-test with precision", pr), { skip_if_no_gpu() for (lag in 1:2) { g <- gpuGranger(x, lag = lag, precision = pr) for (i in 1:3) for (j in 1:3) if (i != j) { ref <- grangerRef(x[, j], x[, i], lag) - if (pr == "double") { - expect_equal(g$fStatistics[i, j], unname(ref["f"]), tolerance = 1e-8) - expect_equal(g$pValues[i, j], unname(ref["p"]), tolerance = 1e-6) - } else { - # near-zero F loses relative precision to float32 cancellation in - # the residual sums of squares; compare on an absolute floor - expect_lt(abs(g$fStatistics[i, j] - ref["f"]), - 1e-3 * max(1, ref["f"])) - expect_lt(abs(g$pValues[i, j] - ref["p"]), 1e-3) - } + expect_equal(g$fStatistics[i, j], unname(ref["f"]), tolerance = 1e-8) + expect_equal(g$pValues[i, j], unname(ref["p"]), tolerance = 1e-8) } } }) - test_that(paste("gpuGranger x/y form agrees with the single-matrix form in", + test_that(paste("gpuGranger x/y form agrees with the single-matrix form,", pr), { skip_if_no_gpu() g <- gpuGranger(x, lag = 2, precision = pr) # several x columns, one y column gxy <- gpuGranger(x[, 2:3], x[, 1, drop = FALSE], lag = 2, precision = pr) expect_equal(gxy$fStatistics[, 1], g$fStatistics[2:3, 1], - tolerance = tolFor(pr), ignore_attr = TRUE) + tolerance = 1e-10, ignore_attr = TRUE) # one x column, several y columns gxy <- gpuGranger(x[, 1, drop = FALSE], x[, 2:3], lag = 2, precision = pr) expect_equal(gxy$fStatistics[1, ], g$fStatistics[1, 2:3], - tolerance = tolFor(pr), ignore_attr = TRUE) + tolerance = 1e-10, ignore_attr = TRUE) }) } +test_that("gpuGranger is accurate on long random walks and trends", { + skip_if_no_gpu() + set.seed(7) + n <- 5000 + walks <- apply(matrix(rnorm(n * 4), n, 4), 2, cumsum) + walks[2:n, 3] <- walks[2:n, 3] + 0.05 * walks[1:(n - 1), 1] + drift <- apply(1 + matrix(rnorm(n * 3), n, 3), 2, cumsum) + trend <- 0.01 * seq_len(n) + matrix(rnorm(n * 3), n, 3) + for (lag in c(1, 3)) for (x in list(walks, drift, trend)) { + g <- expect_no_warning(gpuGranger(x, lag = lag)) + p <- ncol(x) + for (i in 1:p) for (j in 1:p) if (i != j) { + ref <- grangerRef(x[, j], x[, i], lag) + expect_equal(g$fStatistics[i, j], unname(ref["f"]), tolerance = 1e-7, + info = paste("pair", i, j, "lag", lag)) + } + } +}) + +test_that("gpuGranger resolves a cause nearly equal to the effect", { + skip_if_no_gpu() + set.seed(50) + e <- cumsum(rnorm(2000)) + x <- cbind(e, e + 1e-6 * rnorm(2000), cumsum(rnorm(2000))) + g <- expect_no_warning(gpuGranger(x, lag = 2)) + for (ij in list(c(2, 1), c(3, 1))) { + ref <- grangerRef(x[, ij[2]], x[, ij[1]], 2, tol = 1e-12) + expect_equal(g$fStatistics[ij[1], ij[2]], unname(ref["f"]), + tolerance = 1e-6) + } +}) + +test_that("gpuGranger flags what it cannot resolve, and only that", { + skip_if_no_gpu() + set.seed(51) + x <- apply(matrix(rnorm(150 * 3), 150, 3), 2, cumsum) + # a cause that is the effect: nothing of its lags is left beside them + expect_warning(g <- gpuGranger(x, x[, 2, drop = FALSE], lag = 2), + "1 pair\\(s\\) are too ill-conditioned") + expect_true(is.na(g$fStatistics[2, 1])) + expect_true(is.finite(g$fStatistics[1, 1]) && is.finite(g$fStatistics[3, 1])) + # a constant series has no tests, and no warning + xc <- x + xc[, 2] <- 4 + g <- expect_no_warning(gpuGranger(xc, lag = 2)) + expect_true(all(is.na(g$fStatistics[2, ])) && all(is.na(g$fStatistics[, 2]))) + expect_true(is.finite(g$fStatistics[1, 3])) + # series integrated three times at a long lag: every pair is either + # flagged (NA, with the warning) or within 0.1 percent of least squares + set.seed(1) + x3 <- matrix(rnorm(3000 * 3), 3000, 3) + for (k in 1:3) x3 <- apply(x3, 2, cumsum) + warned <- 0 + g <- withCallingHandlers(gpuGranger(x3, lag = 5), + warning = function(w) { + warned <<- warned + 1 + invokeRestart("muffleWarning") + }) + off <- row(g$fStatistics) != col(g$fStatistics) + flagged <- sum(is.na(g$fStatistics[off])) + expect_equal(warned, as.integer(flagged > 0)) + for (i in 1:3) for (j in 1:3) if (i != j && is.finite(g$fStatistics[i, j])) { + ref <- grangerRef(x3[, j], x3[, i], 5, tol = 1e-12) + expect_equal(g$fStatistics[i, j], unname(ref["f"]), tolerance = 1e-3, + info = paste("pair", i, j)) + } +}) + +test_that("gpuGranger agrees with lm() across several tiles of effects", { + skip_if_no_gpu() + # n = 20,000 at lag 3 puts about 34 effects in a tile, so 50 series + # take two tiles: the second tile's indexing, its diagonal and the x/y + # form on a subset are all exercised + set.seed(8) + n <- 20000 + x <- matrix(rnorm(n * 50), n, 50) + x[2:n, 45] <- x[2:n, 45] + 0.1 * x[1:(n - 1), 3] + g <- expect_no_warning(gpuGranger(x, lag = 3)) + expect_true(all(is.na(diag(g$fStatistics)))) + expect_true(all(is.finite(g$fStatistics[row(g$fStatistics) != col(g$fStatistics)]))) + pairs <- rbind(c(3, 45), c(45, 3), c(1, 50), c(50, 1), c(34, 35), c(35, 34), + c(10, 40), c(40, 49), c(49, 48), c(2, 36)) + for (r in seq_len(nrow(pairs))) { + i <- pairs[r, 1]; j <- pairs[r, 2] + ref <- grangerRef(x[, j], x[, i], 3) + expect_equal(g$fStatistics[i, j], unname(ref["f"]), tolerance = 1e-7, + info = paste("pair", i, j)) + } + expect_gt(g$fStatistics[3, 45], 20) + gxy <- gpuGranger(x[, 1:5], x[, 40:50], lag = 3) + expect_equal(gxy$fStatistics, g$fStatistics[1:5, 40:50], tolerance = 1e-10, + ignore_attr = TRUE) +}) + +test_that("gpuGranger checks its arguments", { + skip_if_no_gpu() + expect_error(gpuGranger(x, lag = 0), "positive whole number") + expect_error(gpuGranger(x, lag = 1.5), "positive whole number") + expect_error(gpuGranger(x, lag = 17), "at most 16") + expect_error(gpuGranger(x, lag = 1e10), "at most 16") + expect_error(gpuGranger(x[, 0], lag = 1), "at least one column") + expect_error(gpuGranger(x[1:7, ], lag = 2), "too short") + expect_error(gpuGranger(x, x[1:100, ], lag = 2), "different numbers of rows") + expect_equal(dim(gpuGranger(x, lag = 16)$fStatistics), c(3, 3)) +}) + test_that("gpuGranger's result says which way it reads", { skip_if_no_gpu() g <- gpuGranger(x, lag = 2) diff --git a/vignettes/benchmarks.Rmd b/vignettes/benchmarks.Rmd index ba59b6c..e810aa7 100644 --- a/vignettes/benchmarks.Rmd +++ b/vignettes/benchmarks.Rmd @@ -91,9 +91,14 @@ decomposing large matrices (`gpuSolve`, `gpuPrcomp`). **Single or double precision.** Consumer and workstation GPUs run double precision many times slower than single; data-centre GPUs do not. For -statistics whose result is a correlation, a distance or a test statistic, -single precision usually loses nothing that matters; for least squares on -badly conditioned data, use double. +statistics whose result is a correlation or a distance, single precision +usually loses nothing that matters; for least squares, and for anything +that subtracts two nearly equal sums of squares, use double. `gpuGranger` +is the clear case: its F statistic is a difference of two residual sums +of squares, and in single precision it was off by more than 10 percent on +long random walks, so it now always computes in double, whatever +`precision` says. (The Granger timings in this vignette were measured +before that change; `bench/run-benchmarks.R` repeats them.) ## When the GPU pays off