Добавление скаляра к матрице эффективно в Julia

Question

Мне нужно добавить скаляр ко всем элементам огромной матрицы. Матрица будет как можно больше. В примере я буду использовать размер 2 GiB, но в моем реальном вычислении он будет намного больше.

A = rand(2^14, 2^14) 

Если я исполняю

A += 1.0 

Юлия выделяет дополнительные 2 Гб памяти. Операция занимает около 1 секунды. Я мог бы использовать for цикл:

for jj = 1:size(A, 2), ii = 1:size(A, 1) 
    A[ii, jj] = A[ii, jj] + 1.0 
end 

Это не выделяет память, но она занимает одну минуту. Оба подхода для меня нежизнеспособны, потому что первый нарушает ограничения памяти, а второй явно неэффективен. Для элементарного умножения существует scal!, в котором используется BLAS. Есть ли способ выполнить добавление так же эффективно, как умножение, используя scal!?

Answer 1

@ Ответ DSM является хорошим. Здесь есть несколько вещей, которые я хотел бы затронуть, кроме того. Причина, по которой ваш цикл слишком медленный, заключается в том, что A является непостоянной глобальной переменной, и ваш код напрямую мутирует этот глобальный. Поскольку A непостоянен, код должен защищать от возможности A, становясь другим значением с другим типом в любой момент во время выполнения цикла. Код должен искать тип и расположение A на каждой итерации цикла и динамически отправлять вызовы методов в выражении A[ii, jj] = A[ii, jj] + 1.0 - это вызов getindex, + и setindex!, все из которых зависят от статически неизвестного типа A , Вы можете сразу же получить гораздо более высокую производительность, просто делаю эту работу в функции:

julia> A = rand(2^10, 2^10); 

julia> @time for jj = 1:size(A, 2), ii = 1:size(A, 1) 
      A[ii, jj] += 1 
     end 
elapsed time: 0.288340785 seconds (84048040 bytes allocated, 15.59% gc time) 

julia> function inc!(A) 
      for jj = 1:size(A, 2), ii = 1:size(A, 1) 
       A[ii, jj] += 1 
      end 
     end 
inc! (generic function with 1 method) 

julia> @time inc!(A) 
elapsed time: 0.006076414 seconds (171336 bytes allocated) 

julia> @time inc!(A) 
elapsed time: 0.000888457 seconds (80 bytes allocated) 

избежать неконстантых глобал как это первая рекомендация в Performance Tips разделе руководства. Возможно, вам захочется ознакомиться с остальной частью этой главы.

Далее мы можем улучшить работу функции inc! с помощью @inbounds аннотацию, чтобы указать, что ограничивает проверка не является необходимой для этого кода, а также с использованием линейной индексации вместо двумерный индексации:

julia> function inc!(A) 
      @inbounds for i = 1:length(A) 
       A[i] += 1 
      end 
     end 
inc! (generic function with 1 method) 

julia> @time inc!(A) 
elapsed time: 0.000637934 seconds (80 bytes allocated) 

Большая часть ускорения связана с аннотацией @inbounds, а не с линейной индексацией, хотя это дает небольшое ускорение скорости. Однако аннотацию @inbounds следует использовать экономно, и только там, где каждый уверен, что индексирование не может быть вне пределов и производительности, имеет первостепенное значение. Как вы можете видеть, дополнительное улучшение производительности, хотя и существующее, не является ошеломляющим. Большая часть выгоды исходит не от прямого изменения глобальных переменных.

Answer 2

Вы могли бы сделать на месте операции вещания:

julia> A = rand(2^14, 2^14); A[1:5, 1:5] 
5x5 Array{Float64,2}: 
0.229662 0.680236 0.131202 0.111664 0.802698 
0.500575 0.580994 0.385844 0.983806 0.324382 
0.701694 0.577749 0.532591 0.0508955 0.94325 
0.592929 0.00319653 0.759241 0.448704 0.706204 
0.867945 0.0413606 0.586151 0.82561 0.679233 

julia> @time broadcast!(.+, A, A, 100); 
elapsed time: 0.382669486 seconds (11490976 bytes allocated) 

julia> A[1:5, 1:5] 
5x5 Array{Float64,2}: 
100.23 100.68 100.131 100.112 100.803 
100.501 100.581 100.386 100.984 100.324 
100.702 100.578 100.533 100.051 100.943 
100.593 100.003 100.759 100.449 100.706 
100.868 100.041 100.586 100.826 100.679 

, который использует только в общей сложности ~ 2G памяти.