Ускорение в 14 000 раз или Победа компьютерной науки

Как разработчику научного ПО мне приходится много программировать. И большинство людей из других научных областей склонны думать, что программирование — это «просто» набросать код и запустить его. У меня хорошие рабочие отношения со многими коллегами, в том числе из других стран… Физика, климатология, биология и т. д. Но когда дело доходит до разработки ПО, то складывается отчётливое впечатление, что они думают: «Эй, что тут может быть сложного?! Мы просто записываем несколько инструкций о том, что должен сделать компьютер, нажимаем кнопку „Выполнить” и готово — получаем ответ!»

Проблема в том, что невероятно легко написать инструкции, которые означают не то, что вы думаете. Например, программа может совершенно не поддаваться интерпретации компьютером. Кроме того, нет буквально никакого способа определить, завершится ли программа вообще, не выполнив её. И есть много, очень много способов сильно замедлить выполнение программы. В смысле… реально замедлить. Так замедлить, что выполнение займёт всю вашу жизнь или больше. Это чаще всего происходит с программами, которые написаны людьми без компьютерного образования, то есть учёными из других областей. Моя работа — исправлять такие программы.

Люди не понимают, что информатика учит вас теории вычислений, сложности алгоритмов, вычислимости (то есть можем ли мы действительно что-то вычислить? Слишком часто мы считаем само собой разумеющимся, что можем!) Информатика даёт знания, логику и методы анализа, помогающие написать код, который выполнится за минимальное количество времени или с минимальным использованием ресурсов.

Позвольте показать пример огромной оптимизации одного простого скрипта, написанного моим коллегой.

В климатологии мы много масштабируем. Мы берём показания температуры и осадков из крупномасштабной глобальной климатической модели и сопоставляем их с мелкомасштабной локальной сеткой. Допустим, глобальная сетка равна 50×25, а локальная — 1000×500. Для каждой ячейки сетки в локальной сетке мы хотим знать, какой ячейке сетки в глобальной сетке она соответствует.

Простой способ представить проблему — это минимизировать расстояние между L[n] и G[n]. Получается такой поиск:

для каждой локальной ячейки L[i]:
  для каждой глобальной ячейки G[j]:
     вычислить расстояние между L[i] и G[j]
  найти минимальное расстояние в наборе L[i] * G
  вернуть индекс минимума

Кажется достаточно простым. Однако если посмотреть внимательнее, то мы делаем много дополнительной работы. Посмотрите на алгоритм с точки зрения размера входных данных.

для каждой локальной ячейки L[i]:                           # Сделать L раз
  для каждой глобальной ячейки G[j]:                        # Сделать L x G раз
     вычислить расстояние между L[i] и G[j]                 # Сделать L x G раз
  найти минимальное расстояние в наборе d[i*j]              # Прочитать G ячеек L раз (L x G)
  найти индекс, для которого ячейка соответствует минимуму  # Прочитать G ячеек L раз (L x G)

Код выглядит примерно так:

obs.lon <- ncvar_get(nc.obs, 'lon')
obs.lat <- ncvar_get(nc.obs, 'lat')
n.lon <- length(obs.lon)
n.lat <- length(obs.lat)

obs.lats <- matrix(obs.lat, nrow=n.lon, ncol=n.lat, byrow=TRUE)
obs.lons <- matrix(obs.lon, nrow=n.lon, ncol=n.lat)
obs.time <- netcdf.calendar(nc.obs)

gcm.lon <- ncvar_get(nc.gcm, 'lon')-360
gcm.lat <- ncvar_get(nc.gcm, 'lat')
gcm.lats <- matrix(gcm.lat, ncol=length(gcm.lat), nrow=length(gcm.lon),
                   byrow=TRUE)
gcm.lons <- matrix(gcm.lon, ncol=length(gcm.lat), nrow=length(gcm.lon))
gcm.lons.lats <- cbind(c(gcm.lons), c(gcm.lats))

# Figure out which GCM grid boxes are associated with each fine-scale grid point
# Confine search to 10 deg. x 10 deg. neighbourhood

dxy <- 10
mdist <- function(x, y)
    apply(abs(sweep(data.matrix(y), 2, data.matrix(x), '-')), 1, sum)
nn <- list()
for (i in seq_along(obs.lons)) {
    if((i %% 500)==0) cat(i, '')
    gcm.lims <- ((gcm.lons.lats[,1] >= (obs.lons[i]-dxy)) &
                 (gcm.lons.lats[,1] <= (obs.lons[i]+dxy))) &
                ((gcm.lons.lats[,2] >= (obs.lats[i]-dxy)) &
                 (gcm.lons.lats[,2] <= (obs.lats[i]+dxy)))
    gcm.lims <- which(gcm.lims)
    nn.min <- which.min(mdist(c(obs.lons[i], obs.lats[i]),
                        gcm.lons.lats[gcm.lims,]))
    nn[[i]] <- gcm.lims[nn.min]
}
nn <- unlist(nn)

Похоже на простой алгоритм. «Просто» вычислить расстояния, а затем найти минимум. Но в такой реализации по мере роста числа локальных ячеек наша стоимость вычислений растёт на её произведение с числом ячеек глобальной сетки. Канадские данные ANUSPLIN содержат 1068×510 ячеек (в общей сложности 544 680). Предположим, что наш GCM содержит 50×25 ячеек (в общей сложности 1250). Таким образом, стоимость внутреннего цикла в «каких-то вычислительных единицах» равна:

где члены — это константы, соответствующие стоимости вычисления расстояния между двумя точками, нахождения минимальной точки и нахождения индекса массива. На самом деле, мы не заботимся (особо) о постоянных членах, потому что они не зависят от размера входных данных. Так что можно просто сложить их и оценить стоимость вычисления;

Таким образом, для этого набора входных данных наша стоимость составляет .

680 миллионов.

Кажется, много… хотя, кто его знает? Компьютеры ведь быстрые, верно? Если мы запустим наивную реализацию, в реальности она отработает чуть быстрее 1668 секунд, что составляет чуть меньше получаса.

> source('BCCA/naive.implementation.R')
500 1000 1500 2000 2500 3000 ... 543000 543500 544000 544500 [1] "Elapsed Time"
    user   system  elapsed 
1668.868    8.926 1681.728 

Но должна ли программа выполняться 30 минут? Вот в чём дело. Мы сравниваем две сетки, в каждой из которых есть масса структур, которые мы никак не задействовали. Например, широты и долготы в обеих сетках расположены в отсортированном порядке. Поэтому, чтобы найти число, не нужно просматривать весь массив. Вы можете использовать алгоритм деления пополам — смотрите точку в середине и решаете, какая половина массива нам нужна. Тогда поиск по всему пространству займёт логарифм по основанию два от всего пространства поиска.

Другая важная структура, которой мы не воспользовались — тот факт, что широты идут в измерении , а долготы — в измерении . Таким образом, вместо выполнения операции раз мы можем выполнить её раз. Это огромная оптимизация.

Как это выглядит в псевдокоде?

Для каждого local[x]:
    bisect_search(local[x], Global[x])

Для каждого local[y]:
    bisect_search(local[y], Global[y])

возвращается 2d сетка результатов поиска для каждого измерения

В коде:

## Perform a binary search on the *sorted* vector v
## Return the array index of the element closest to x
find.nearest <- function(x, v) {
    if (length(v) == 1) {
        return(1)
    }
    if (length(v) == 2) {
        return(which.min(abs(v - x)))
    }
    mid <- ceiling(length(v) / 2)
    if (x == v[mid]) {
        return(mid)
    } else if (x < v[mid]) {
        return(find.nearest(x, v[1:mid]))
    }
    else {
        return((mid - 1) + find.nearest(x, v[mid:length(v)]))
    }
}

regrid.one.dim <- function(coarse.points, fine.points) {
    return(sapply(fine.points, find.nearest, coarse.points))
}

## Take a fine scale (e.g. ANUSPLINE) grid of latitudes and longitudes
## and find the indicies that correspond to a coarse scale (e.g. a GCM) grid
## Since the search is essentially a minimizing distance in 2 dimensions
## We can actually search independently in each dimensions separately (which
## is a huge optimization, making the run time x + y instead of x * y) and
## then reconstruct the indices to create a full grid
regrid.coarse.to.fine <- function(coarse.lats, coarse.lons, fine.lats, fine.lons) {
    xi <- regrid.one.dim(gcm.lon, obs.lon)
    yi <- regrid.one.dim(gcm.lat, obs.lat)
    ## Two dimensional grid of indices
    xi <- matrix(xi, ncol=length(fine.lats), nrow=length(fine.lons), byrow=F)
    yi <- matrix(yi, ncol=length(fine.lats), nrow=length(fine.lons), byrow=T)
    return(list(xi=xi, yi=yi))
}

Стоимость каждого поиска методом деления пополам — это log от входного размера. На этот раз наш входной размер разделён на пространство X и Y, поэтому будем использовать и для Global, Local, X и Y.

Затраты получаются в 553 076. Что ж, 553 тысячи звучит намного лучше, чем 680 миллионов. Как это отразится на времени выполнения?

> ptm <- proc.time(); rv <- regrid.coarse.to.fine(gcm.lat, gcm.lon, obs.lat, obs.lon); print('Elapsed Time'); print(proc.time() - ptm)[1] "Elapsed Time"
   user  system elapsed 
  0.117   0.000   0.117 
> str(rv)
List of 2
 $ xi: num [1:1068, 1:510] 15 15 15 15 15 15 15 15 15 15 ...
 $ yi: num [1:1068, 1:510] 13 13 13 13 13 13 13 13 13 13 ...
> 

0,117 секунды. То, что раньше выполнялось почти полчаса, теперь занимает чуть больше одной десятой секунды.

> 1668.868 / .117
[1] 14263.83

Таааааак… Понимаю, что я обучен заниматься такой оптимизацией и знаю, как её делать хорошо, но даже меня удивляет и впечатляет, насколько значительно ускорение. Примерно в 14 тысяч раз.

Раньше скрипт работал настолько медленно, что сохранял результат на диск для проверки учёным вручную, прежде чем продолжить. Теперь всё вычисляется в мгновение ока. Такие расчёты мы проводим сотни раз, так что в итоге экономим дни и даже недели вычислительного времени. И получаем возможность лучше взаимодействовать с системой, так что рабочее время учёных проходит с большей пользой… они не сидят без дела, ожидая окончания вычислений. Всё готово сразу.

Должен подчеркнуть, что для этих эпических улучшений производительности не требуется покупка каких-то больших компьютерных систем, распараллеливание или увеличение сложности… на самом деле код более быстрого алгоритма даже более простой и универсальный! Эта полная победа достигнута просто благодаря вдумчивому чтению кода и наличию некоторых знаний о вычислительной сложности.