Рекомендации по производительности

В следующих разделах мы кратко рассмотрим несколько техник, которые могут помочь сделать ваш код на Julia максимально быстрым.

Критически важный для производительности код должен находиться внутри функции

Любой код, критически важный для производительности, должен находиться внутри функции. Код внутри функций, как правило, выполняется намного быстрее, чем код верхнего уровня, из-за того, как работает компилятор Julia.

Использование функций важно не только для производительности: функции более переиспользуемы и проверяемы, и они поясняют, какие шаги выполняются и каковы их входные и выходные данные. Написание функций, а не просто скриптов — это также рекомендация руководства по стилю Julia.

Функции должны принимать аргументы, а не работать непосредственно с глобальными переменными, см. следующий пункт.

Избегайте нетипизированных глобальных переменных

Значение нетипизированной глобальной переменной может измениться в любой момент, что, возможно, приведёт к изменению её типа. Это затрудняет компилятору оптимизировать код, использующий глобальные переменные. Это также относится к переменным типа, т. е. алиасам типов на глобальном уровне. Переменные должны быть локальными или передаваться в качестве аргументов функциям, когда это возможно.

Мы обнаружили, что глобальные имена часто являются константами, и объявление их как таковых значительно улучшает производительность:

const DEFAULT_VAL = 0

Если известно, что глобальная переменная всегда имеет один и тот же тип, тип должен быть аннотирован.

Использование нетипизированных глобальных переменных можно оптимизировать, аннотируя их типы в месте использования:

global x = rand(1000)

function loop_over_global()
    s = 0.0
    for i in x::Vector{Float64}
        s += i
    end
    return s
end

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

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

В следующем сеансе REPL:

julia> x = 1.0

эквивалентно:

julia> global x = 1.0

поэтому все проблемы производительности, обсуждавшиеся ранее, применяются.

Измеряйте производительность с помощью `@time` и обращайте внимание на выделение памяти

Полезным инструментом для измерения производительности является макрос @time. Мы здесь повторяем пример с глобальной переменной выше, но на этот раз с удалённой аннотацией типа:

julia> x = rand(1000);

julia> function sum_global()
           s = 0.0
           for i in x
               s += i
           end
           return s
       end;

julia> @time sum_global()
  0.011539 seconds (9.08 k allocations: 373.386 KiB, 98.69% compilation time)
523.0007221951678

julia> @time sum_global()
  0.000091 seconds (3.49 k allocations: 70.156 KiB)
523.0007221951678

При первом вызове (@time sum_global()) функция компилируется. (Если вы ещё не использовали @time в этом сеансе, это также скомпилирует функции, необходимые для измерения времени.) Вы не должны принимать результаты этого запуска всерьёз. Для второго запуска обратите внимание, что помимо отчёта о времени, также указано, что было выделено значительное количество памяти. Здесь мы просто вычисляем сумму всех элементов в векторе 64-битных чисел с плавающей точкой, поэтому выделение (кучевой) памяти не должно потребоваться.

Мы должны уточнить, что то, что @time сообщает, — это конкретно выделение кучевой памяти, которое обычно необходимо для мутабельных объектов или для создания/увеличения контейнеров переменного размера (таких как Array или Dict, строки или объекты с «нестабильным типом», тип которых известен только во время выполнения). Выделение (или освобождение) таких блоков памяти может потребовать дорогостоящего вызова libc (например, через malloc на C), и их необходимо отслеживать для сбора мусора. В отличие от этого, неизменяемые значения, такие как числа (кроме больших чисел), кортежи и неизменяемые struct можно хранить гораздо дешевле, например, в стеке или памяти процессора, поэтому обычно не нужно беспокоиться о затратах на производительность при их «выделении».

Неожиданное выделение памяти почти всегда является признаком какой-то проблемы с вашим кодом, обычно проблемы со стабильностью типа или создания многих небольших временных массивов. Следовательно, помимо самого выделения, очень вероятно, что сгенерированный для вашей функции код далёк от оптимального. Воспринимайте такие указания всерьёз и следуйте приведённым ниже рекомендациям.

В этом конкретном случае выделение памяти вызвано использованием глобальной переменной с нестабильным типом x, поэтому, если вместо этого передать x в качестве аргумента функции, она больше не выделяет память (оставшееся выделение, указанное ниже, вызвано выполнением макроса @time в глобальной области видимости) и работает значительно быстрее после первого вызова:

julia> x = rand(1000);

julia> function sum_arg(x)
           s = 0.0
           for i in x
               s += i
           end
           return s
       end;

julia> @time sum_arg(x)
  0.007551 seconds (3.98 k allocations: 200.548 KiB, 99.77% compilation time)
523.0007221951678

julia> @time sum_arg(x)
  0.000006 seconds (1 allocation: 16 bytes)
523.0007221951678

Единственное выделение, видимое, — это выделение, вызванное самим макросом @time в глобальной области видимости. Если мы вместо этого запустим измерение времени в функции, мы увидим, что выделение действительно не выполняется:

julia> time_sum(x) = @time sum_arg(x);

julia> time_sum(x)
  0.000002 seconds
523.0007221951678

В некоторых ситуациях вашей функции может потребоваться выделение памяти в процессе её работы, и это может усложнить простую картину, описанную выше. В таких случаях рассмотрите возможность использования одного из инструментов ниже для диагностики проблем или написания версии вашей функции, которая отделяет выделение от её алгоритмических аспектов (см. Предварительное выделение выходов).

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

Инструменты

Julia и её экосистема пакетов включают инструменты, которые могут помочь вам диагностировать проблемы и улучшить производительность вашего кода:

Профилирование позволяет измерить производительность вашего работающего кода и выявить строки, которые являются узкими местами. Для сложных проектов пакет ProfileView может помочь вам визуализировать результаты профилирования.
Пакет JET может помочь вам найти распространённые проблемы производительности в вашем коде.
Неожиданно большие выделения памяти — как сообщается @time, @allocated или профилировщиком (через вызовы процедур сбора мусора) — указывают на то, что могут быть проблемы с вашим кодом. Если вы не видите другой причины для выделений, подозревайте проблему с типом. Вы также можете запустить Julia с опцией --track-allocation=user и просмотреть созданные файлы *.mem для получения информации о том, где происходят эти выделения. См. Анализ выделения памяти.
@code_warntype генерирует представление вашего кода, которое может быть полезно для поиска выражений, приводящих к неопределённости типа. См. @code_warntype ниже.

Избегайте контейнеров с параметрами абстрактного типа

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

Рассмотрим следующее:

julia> a = Real[]
Real[]

julia> push!(a, 1); push!(a, 2.0); push!(a, π)
3-element Vector{Real}:
 1
 2.0
 π = 3.1415926535897...

Поскольку a является массивом абстрактного типа Real, он должен быть способен содержать любое значение Real. Поскольку объекты Real могут иметь произвольный размер и структуру, a должен быть представлен в виде массива указателей на индивидуально выделенные объекты Real. Однако, если мы вместо этого разрешим хранить только числа одного типа, например, Float64, в a они могут храниться более эффективно:

julia> a = Float64[]
Float64[]

julia> push!(a, 1); push!(a, 2.0); push!(a,  π)
3-element Vector{Float64}:
 1.0
 2.0
 3.141592653589793

Присвоение чисел в a теперь преобразует их в Float64, а a будет храниться как непрерывный блок 64-битных значений с плавающей точкой, которые могут быть эффективно обработаны.

Если вы не можете избежать контейнеров с абстрактными типами значений, иногда лучше параметризовать с Any, чтобы избежать проверки типа во время выполнения. Например, IdDict{Any, Any} работает лучше, чем IdDict{Type, Vector}

См. также обсуждение в разделе Параметризованные типы.

Объявления типов

Во многих языках с необязательными объявлениями типов добавление объявлений — основной способ ускорить работу кода. Это не относится к Julia. В Julia компилятор, как правило, знает типы всех аргументов функций, локальных переменных и выражений. Однако есть несколько конкретных случаев, когда объявления полезны.

Избегайте полей с абстрактным типом

Типы могут быть объявлены без указания типов их полей:

julia> struct MyAmbiguousType
           a
       end

Это позволяет a быть любого типа. Это часто бывает полезно, но у этого есть недостаток: для объектов типа MyAmbiguousType, компилятор не сможет сгенерировать код высокой производительности. Причина в том, что компилятор использует типы объектов, а не их значения, чтобы определить, как построить код. К сожалению, о объекте типа MyAmbiguousType можно вывести очень мало:

julia> b = MyAmbiguousType("Hello")
MyAmbiguousType("Hello")

julia> c = MyAmbiguousType(17)
MyAmbiguousType(17)

julia> typeof(b)
MyAmbiguousType

julia> typeof(c)
MyAmbiguousType

Значения b и c имеют одинаковый тип, но их внутреннее представление данных в памяти сильно различается. Даже если вы хранили только числовые значения в поле a, тот факт, что представление в памяти UInt8 отличается от Float64, также означает, что процессору необходимо обрабатывать их с помощью двух разных типов инструкций. Поскольку необходимая информация недоступна в типе, такие решения должны приниматься во время выполнения. Это снижает производительность.

Вы можете улучшить ситуацию, объявив тип a. Здесь мы сосредоточены на случае, когда a может быть любым из нескольких типов, в этом случае естественным решением является использование параметров. Например:

julia> mutable struct MyType{T<:AbstractFloat}
           a::T
       end

Это лучший выбор, чем

julia> mutable struct MyStillAmbiguousType
           a::AbstractFloat
       end

потому что в первом варианте тип a определяется по типу объекта-оболочки. Например:

julia> m = MyType(3.2)
MyType{Float64}(3.2)

julia> t = MyStillAmbiguousType(3.2)
MyStillAmbiguousType(3.2)

julia> typeof(m)
MyType{Float64}

julia> typeof(t)
MyStillAmbiguousType

Тип поля a легко определяется по типу m, но не по типу t. Действительно, в t можно изменить тип поля a:

julia> typeof(t.a)
Float64

julia> t.a = 4.5f0
4.5f0

julia> typeof(t.a)
Float32

В отличие от этого, после создания m, тип m.a изменить нельзя:

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float64

Тот факт, что тип m.a известен по типу m, а также то, что его тип не может измениться внутри функции, позволяет компилятору генерировать высокооптимизированный код для объектов, таких как m, но не для объектов, таких как t.

Конечно, всё это верно только в том случае, если мы создаём m с конкретным типом. Мы можем нарушить это, явно создав его с абстрактным типом:

julia> m = MyType{AbstractFloat}(3.2)
MyType{AbstractFloat}(3.2)

julia> typeof(m.a)
Float64

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float32

Для всех практических целей такие объекты ведут себя идентично объектам MyStillAmbiguousType.

Очень показательно сравнить объём кода, генерируемого для простой функции

func(m::MyType) = m.a+1

используя

code_llvm(func, Tuple{MyType{Float64}})
code_llvm(func, Tuple{MyType{AbstractFloat}})

По причинам объёма результаты здесь не показаны, но вы можете попробовать сделать это самостоятельно. Поскольку тип полностью определён в первом случае, компилятору не нужно генерировать код для разрешения типа во время выполнения. Это приводит к более короткому и быстрому коду.

Также следует помнить, что неполностью параметризованные типы ведут себя как абстрактные типы. Например, хотя полностью определённый Array{T,n} является конкретным, Array без указанных параметров не является конкретным:

julia> !isconcretetype(Array), !isabstracttype(Array), isstructtype(Array), !isconcretetype(Array{Int}), isconcretetype(Array{Int,1})
(true, true, true, true, true)

В этом случае лучше избегать объявления MyType с полем a::Array, и вместо этого объявить поле как a::Array{T,N} или как a::A, где {T,N} или A являются параметрами MyType.

Предыдущий совет особенно полезен, когда поля структуры предназначены для функций или, более общим образом, для вызываемых объектов. Очень заманчиво определить структуру следующим образом:

struct MyCallableWrapper
    f::Function
end

Но так как Function — абстрактный тип, каждый вызов wrapper.f потребует динамической диспатчизации из-за неустойчивости типа при доступе к полю f. Вместо этого вы должны написать что-то вроде:

struct MyCallableWrapper{F}
    f::F
end

что имеет почти идентичное поведение, но будет значительно быстрее (потому что устранена неустойчивость типа). Обратите внимание, что мы не накладываем F<:Function: это означает, что вызываемые объекты, не являющиеся подтипами Function, также разрешены для поля f.

Избегайте полей с абстрактными контейнерами

Такие же рекомендации действительны и для типов контейнеров:

julia> struct MySimpleContainer{A<:AbstractVector}
           a::A
       end

julia> struct MyAmbiguousContainer{T}
           a::AbstractVector{T}
       end

julia> struct MyAlsoAmbiguousContainer
           a::Array
       end

Например:

julia> c = MySimpleContainer(1:3);

julia> typeof(c)
MySimpleContainer{UnitRange{Int64}}

julia> c = MySimpleContainer([1:3;]);

julia> typeof(c)
MySimpleContainer{Vector{Int64}}

julia> b = MyAmbiguousContainer(1:3);

julia> typeof(b)
MyAmbiguousContainer{Int64}

julia> b = MyAmbiguousContainer([1:3;]);

julia> typeof(b)
MyAmbiguousContainer{Int64}

julia> d = MyAlsoAmbiguousContainer(1:3);

julia> typeof(d), typeof(d.a)
(MyAlsoAmbiguousContainer, Vector{Int64})

julia> d = MyAlsoAmbiguousContainer(1:1.0:3);

julia> typeof(d), typeof(d.a)
(MyAlsoAmbiguousContainer, Vector{Float64})

Для MySimpleContainer, объект полностью определён своим типом и параметрами, поэтому компилятор может генерировать оптимизированные функции. В большинстве случаев этого, вероятно, будет достаточно.

Хотя компилятор теперь может выполнять свою работу безупречно, есть случаи, когда вы можете захотеть, чтобы ваш код мог выполнять разные действия в зависимости от типа элемента a. Обычно лучший способ добиться этого — обернуть вашу специфическую операцию (здесь, foo) в отдельную функцию:

julia> function sumfoo(c::MySimpleContainer)
           s = 0
           for x in c.a
               s += foo(x)
           end
           s
       end
sumfoo (generic function with 1 method)

julia> foo(x::Integer) = x
foo (generic function with 1 method)

julia> foo(x::AbstractFloat) = round(x)
foo (generic function with 2 methods)

Это сохраняет простоту, позволяя компилятору генерировать оптимизированный код во всех случаях.

Однако есть случаи, когда вам может потребоваться объявить различные версии внешней функции для разных типов элементов или типов AbstractVector поля a в MySimpleContainer. Вы можете сделать это так:

julia> function myfunc(c::MySimpleContainer{<:AbstractArray{<:Integer}})
           return c.a[1]+1
       end
myfunc (generic function with 1 method)

julia> function myfunc(c::MySimpleContainer{<:AbstractArray{<:AbstractFloat}})
           return c.a[1]+2
       end
myfunc (generic function with 2 methods)

julia> function myfunc(c::MySimpleContainer{Vector{T}}) where T <: Integer
           return c.a[1]+3
       end
myfunc (generic function with 3 methods)

julia> myfunc(MySimpleContainer(1:3))
2

julia> myfunc(MySimpleContainer(1.0:3))
3.0

julia> myfunc(MySimpleContainer([1:3;]))
4

Аннотирование значений, взятых из нетипизированных местоположений

Часто удобно работать со структурами данных, которые могут содержать значения любого типа (массивы типа Array{Any}). Но если вы используете одну из этих структур и знаете тип элемента, это помогает поделиться этими знаниями с компилятором:

function foo(a::Array{Any,1})
    x = a[1]::Int32
    b = x+1
    ...
end

Здесь мы знали, что первый элемент a будет Int32. Такая аннотация имеет дополнительное преимущество: она вызывает ошибку во время выполнения, если значение не имеет ожидаемого типа, потенциально выявляя определённые ошибки раньше.

В случае, если тип a[1] не известен точно, x может быть объявлен через x = convert(Int32, a[1])::Int32. Использование функции convert позволяет a[1] быть любым объектом, преобразуемым в Int32 (например, UInt8 ), тем самым повышая общность кода путём ослабления требования к типу. Обратите внимание, что convert в этом контексте также нуждается в аннотации типа для достижения стабильности типа. Это связано с тем, что компилятор не может вывести тип возвращаемого значения функции, даже convert, если не известны типы всех аргументов функции.

Аннотация типа не улучшит (и может даже помешать) производительности, если тип абстрактный или создаётся во время выполнения. Это связано с тем, что компилятор не может использовать аннотацию для специализации последующего кода, а сама проверка типа занимает время. Например, в коде:

function nr(a, prec)
    ctype = prec == 32 ? Float32 : Float64
    b = Complex{ctype}(a)
    c = (b + 1.0f0)::Complex{ctype}
    abs(c)
end

аннотация c ухудшает производительность. Чтобы написать производительный код, связанный с типами, созданными во время выполнения, используйте технику «барьера функций», описанную ниже, и убедитесь, что созданный тип появляется среди типов аргументов ядра функции, чтобы ядро операций правильно специализировалось компилятором. Например, в фрагменте кода как только b создан, он может быть передан в другую функцию k, ядро. Если, например, функция k объявляет b в качестве аргумента типа Complex{T}, где T — параметр типа, то аннотация типа, появляющаяся в операторе присваивания в k, вида:

c = (b + 1.0f0)::Complex{T}

не ухудшает производительность (но и не помогает), так как компилятор может определить тип c в момент компиляции k.

Обратите внимание на случаи, когда Julia избегает специализации

Как эвристика, Julia избегает автоматической специализации по параметрам типа аргумента в трёх конкретных случаях: Type, Function, и Vararg. Julia всегда специализируется, когда аргумент используется внутри метода, но не если аргумент просто передаётся другой функции. Это обычно не влияет на производительность во время выполнения и улучшает производительность компилятора. Если вы обнаружите, что это влияет на производительность во время выполнения в вашем случае, вы можете запустить специализацию, добавив параметр типа в объявление метода. Вот некоторые примеры:

Это не будет специализироваться:

function f_type(t)  # or t::Type
    x = ones(t, 10)
    return sum(map(sin, x))
end

но это будет:

function g_type(t::Type{T}) where T
    x = ones(T, 10)
    return sum(map(sin, x))
end

Эти не будут специализироваться:

f_func(f, num) = ntuple(f, div(num, 2))
g_func(g::Function, num) = ntuple(g, div(num, 2))

но это будет:

h_func(h::H, num) where {H} = ntuple(h, div(num, 2))

Это не будет специализироваться:

f_vararg(x::Int...) = tuple(x...)

но это будет:

g_vararg(x::Vararg{Int, N}) where {N} = tuple(x...)

Достаточно ввести один параметр типа, чтобы заставить специализацию, даже если другие типы не ограничены. Например, это также будет специализироваться, и это полезно, когда аргументы не все одного типа:

h_vararg(x::Vararg{Any, N}) where {N} = tuple(x...)

Обратите внимание, что @code_typed и аналогичные функции всегда показывают вам специализированный код, даже если Julia обычно не специализирует этот вызов метода. Вам нужно проверить внутренности метода, если вы хотите увидеть, генерируются ли специализации при изменении типов аргументов, т.е., если Base.specializations(@which f(...)) содержит специализации для рассматриваемого аргумента.

Разбиение функций на несколько определений

Написание функции как многих небольших определений позволяет компилятору напрямую вызывать наиболее подходящий код или даже включать его.

Вот пример «сложной функции», которую следует записать как несколько определений:

using LinearAlgebra

function mynorm(A)
    if isa(A, Vector)
        return sqrt(real(dot(A,A)))
    elseif isa(A, Matrix)
        return maximum(svdvals(A))
    else
        error("mynorm: invalid argument")
    end
end

Это можно записать более лаконично и эффективно как:

mynorm(x::Vector) = sqrt(real(dot(x, x)))
mynorm(A::Matrix) = maximum(svdvals(A))

Однако следует отметить, что компилятор довольно эффективно оптимизирует удалённые ветви в коде, написанном как пример mynorm.

Написание «стабильных по типу» функций

Когда это возможно, помогает гарантировать, что функция всегда возвращает значение одного и того же типа. Рассмотрим следующее определение:

pos(x) = x < 0 ? 0 : x

Хотя это кажется достаточно безобидным, проблема в том, что 0 — это целое число (типа Int), а x может быть любого типа. Таким образом, в зависимости от значения x, эта функция может вернуть значение любого из двух типов. Это поведение разрешено и может быть желательным в некоторых случаях. Но это легко исправить следующим образом:

pos(x) = x < 0 ? zero(x) : x

Также существует функция oneunit и более общая функция oftype(x, y), которая возвращает y, преобразованное к типу x.

Избегайте изменения типа переменной

Аналогичная проблема «стабильности типа» существует для переменных, используемых многократно в функции:

function foo()
    x = 1
    for i = 1:10
        x /= rand()
    end
    return x
end

Локальная переменная x начинается как целое число и после одной итерации цикла становится числом с плавающей запятой (результат оператора /). Это затрудняет оптимизацию компилятором тела цикла. Существует несколько возможных исправлений:

Инициализировать x с x = 1.0
Явно указать тип x как x::Float64 = 1
Использовать явное преобразование x = oneunit(Float64)
Инициализировать с первой итерации цикла, чтобы x = 1 / rand(), затем выполнить цикл for i = 2:10

Отдельные ядра функций (иначе, барьеры функций)

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

julia> function strange_twos(n)
           a = Vector{rand(Bool) ? Int64 : Float64}(undef, n)
           for i = 1:n
               a[i] = 2
           end
           return a
       end;

julia> strange_twos(3)
3-element Vector{Int64}:
 2
 2
 2

Это следует написать так:

julia> function fill_twos!(a)
           for i = eachindex(a)
               a[i] = 2
           end
       end;

julia> function strange_twos(n)
           a = Vector{rand(Bool) ? Int64 : Float64}(undef, n)
           fill_twos!(a)
           return a
       end;

julia> strange_twos(3)
3-element Vector{Int64}:
 2
 2
 2

Компилятор Julia специализирует код для типов аргументов на границах функций, поэтому в исходной реализации он не знает тип a во время цикла (поскольку он выбирается случайным образом). Поэтому второй вариант, как правило, быстрее, поскольку внутренний цикл может быть перекомпилирован как часть fill_twos! для различных типов a.

Второй вариант также часто является более стилистичным и может привести к большему повторному использованию кода.

Этот шаблон используется в нескольких местах в Julia Base. Например, см. vcat и hcat в abstractarray.jl или функцию fill!, которую мы могли бы использовать вместо написания собственной функции fill_twos!.

Функции, такие как strange_twos , возникают при работе с данными неопределенного типа, например, данными, загруженными из входного файла, который может содержать целые числа, числа с плавающей точкой, строки или что-то еще.

Типы со значениями в качестве параметров

Предположим, вы хотите создать N-мерный массив, размер которого равен 3 вдоль каждой оси. Такие массивы можно создать следующим образом:

julia> A = fill(5.0, (3, 3))
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

Этот подход работает очень хорошо: компилятор может определить, что A является Array{Float64,2}, поскольку он знает тип значения заполнения (5.0::Float64) и размерность ((3, 3)::NTuple{2,Int}). Это означает, что компилятор может сгенерировать очень эффективный код для любого последующего использования A в той же функции.

Но предположим, что вы хотите написать функцию, которая создает массив 3×3×… в произвольных измерениях; вы можете быть склонны написать функцию

julia> function array3(fillval, N)
           fill(fillval, ntuple(d->3, N))
       end
array3 (generic function with 1 method)

julia> array3(5.0, 2)
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

Это работает, но (как вы можете проверить сами, используя @code_warntype array3(5.0, 2)) проблема в том, что тип результата нельзя вывести: аргумент N является значением типа Int, и вывод типов не (и не может) предсказать его значение заранее. Это означает, что код, использующий результат этой функции, должен быть консервативным, проверяя тип при каждом обращении к A; такой код будет очень медленным.

Сейчас один из очень хороших способов решения таких проблем заключается в использовании техники «барьера функций» . Однако в некоторых случаях вы можете захотеть полностью устранить неустойчивость типа. В таких случаях одним из подходов является передача размерности в качестве параметра, например, через Val{T}() (см. "Типы значений"):

julia> function array3(fillval, ::Val{N}) where N
           fill(fillval, ntuple(d->3, Val(N)))
       end
array3 (generic function with 1 method)

julia> array3(5.0, Val(2))
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

Julia имеет специализированную версию ntuple , которая принимает экземпляр Val{::Int} в качестве второго параметра; передавая N в качестве параметра типа, вы делаете его «значение» известным компилятору. В результате эта версия array3 позволяет компилятору предсказать тип возвращаемого значения.

Однако использование таких техник может быть неожиданно тонким. Например, это не помогло бы, если бы вы вызывали array3 из функции такого вида:

function call_array3(fillval, n)
    A = array3(fillval, Val(n))
end

Здесь вы создали ту же проблему заново: компилятор не может угадать, что такое n, поэтому он не знает тип Val(n). Попытка использовать Val, но неправильно, может легко ухудшить производительность во многих ситуациях. (Только в ситуациях, когда вы эффективно объединяете Val с техникой барьера функций, чтобы сделать ядро функции более эффективным, код, подобный приведенному выше, следует использовать.)

Пример правильного использования Val был бы:

function filter3(A::AbstractArray{T,N}) where {T,N}
    kernel = array3(1, Val(N))
    filter(A, kernel)
end

В этом примере N передается в качестве параметра, поэтому его «значение» известно компилятору. В сущности, Val(T) работает только тогда, когда T либо жёстко закодирован/литерален (Val(3)) или уже указан в области типов.

Опасности злоупотребления множественным диспетчерированием (т.е., больше о типах со значениями в качестве параметров)

Как только человек научится ценить множественное диспетчерирование, возникает понятная тенденция переусердствовать и попытаться использовать его для всего. Например, вы можете представить себе его использование для хранения информации, например,

struct Car{Make, Model}
    year::Int
    ...more fields...
end

а затем диспетчерировать по объектам, таким как Car{:Honda,:Accord}(year, args...).

Это может быть целесообразно, если выполняются следующие условия:

Вам требуется ресурсоёмкая обработка каждого Car, и она становится значительно эффективнее, если вы знаете Make и Model в момент компиляции, и общее количество различных Make или Model , которые будут использоваться, не слишком велико.
У вас есть однородные списки одного и того же типа Car для обработки, поэтому вы можете сохранить их все в Array{Car{:Honda,:Accord},N}.

Когда последнее верно, функция, обрабатывающая такой однородный массив, может быть продуктивно специализирована: Julia знает тип каждого элемента заранее (все объекты в контейнере имеют один и тот же конкретный тип), поэтому Julia может «посмотреть» правильные вызовы методов при компиляции функции (устраняя необходимость проверки во время выполнения) и тем самым генерировать эффективный код для обработки всего списка.

Когда это не так, скорее всего, вы не получите никакой выгоды; хуже того, получившийся «комбинаторный взрыв типов» будет контрпродуктивным. Если items[i+1] имеет другой тип, чем item[i], Julia должна искать тип во время выполнения, искать соответствующий метод в таблицах методов, решать (с помощью пересечения типов), какой из них соответствует, определять, был ли он уже скомпилирован JIT (и делать это, если нет), и затем совершать вызов. По сути, вы просите полную систему типов и механизм JIT-компиляции выполнить примерно то же, что и оператор switch или поиск в словаре в вашем собственном коде.

Некоторые эталонные тесты времени выполнения, сравнивающие (1) диспетчерирование типов, (2) поиск в словаре и (3) оператор «switch», можно найти на форуме рассылки.

Возможно, еще хуже, чем влияние на время выполнения, это влияние на время компиляции: Julia будет компилировать специализированные функции для каждого различного Car{Make, Model}; если у вас есть сотни или тысячи таких типов, то каждая функция, принимающая такой объект в качестве параметра (от пользовательской функции get_year , которую вы можете написать сами, до универсальной функции push! в Julia Base), будет иметь сотни или тысячи вариантов компиляции для неё. Каждый из них увеличивает размер кэша скомпилированного кода, длину внутренних списков методов и т.д. Чрезмерный энтузиазм по поводу значений в качестве параметров может легко привести к растрате огромных ресурсов.

Доступ к массивам в порядке памяти по столбцам

Многомерные массивы в Julia хранятся в порядке следования столбцов. Это означает, что массивы складываются один столбец за другим. Это можно проверить, используя функцию vec или синтаксис [:], как показано ниже (обратите внимание, что массив упорядочен [1 3 2 4], а не [1 2 3 4]):

julia> x = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> x[:]
4-element Vector{Int64}:
 1
 3
 2
 4

Эта соглашение для упорядочивания массивов широко распространено во многих языках, таких как Fortran, Matlab и R (и не только). Альтернативой порядку столбцов является порядок строк, принятый в C и Python (numpy ) среди других языков. Запоминание порядка массивов может иметь существенное влияние на производительность при циклическом переборе массивов. Правило для запоминания состоит в том, что для массивов в порядке столбцов индекс первого столбца меняется быстрее всего. По сути, это означает, что циклический перебор будет быстрее, если внутренний индекс цикла — это первый индекс в выражении среза. Имейте в виду, что индексирование массива с помощью : является неявным циклом, который последовательно обращается ко всем элементам в определенном измерении; например, извлечение столбцов может быть быстрее, чем извлечение строк.

Рассмотрим следующий условный пример. Представьте, что нам нужно написать функцию, которая принимает вектор Vector и возвращает квадратную матрицу Matrix, где либо строки, либо столбцы заполнены копиями входного вектора. Предположим, что порядок заполнения (строки или столбцы) не важен (возможно, остальную часть кода можно легко адаптировать соответственно). Мы можем сделать это по крайней мере четырьмя способами (в дополнение к рекомендуемому вызову встроенной функции repeat):

function copy_cols(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[:, i] = x
    end
    return out
end

function copy_rows(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[i, :] = x
    end
    return out
end

function copy_col_row(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for col = inds, row = inds
        out[row, col] = x[row]
    end
    return out
end

function copy_row_col(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for row = inds, col = inds
        out[row, col] = x[col]
    end
    return out
end

Теперь мы измерим время выполнения каждой из этих функций с использованием того же случайного 10000 на 1 входного вектора:

julia> x = randn(10000);

julia> fmt(f) = println(rpad(string(f)*": ", 14, ' '), @elapsed f(x))

julia> map(fmt, [copy_cols, copy_rows, copy_col_row, copy_row_col]);
copy_cols:    0.331706323
copy_rows:    1.799009911
copy_col_row: 0.415630047
copy_row_col: 1.721531501

Обратите внимание, что copy_cols намного быстрее, чем copy_rows. Это ожидаемо, потому что copy_cols учитывает расположение памяти массива по столбцам и заполняет его один столбец за другим. Кроме того, copy_col_row намного быстрее, чем copy_row_col, потому что он следует нашему правилу, согласно которому первый элемент, появляющийся в выражении среза, должен быть связан с внутренним циклом.

Предварительное выделение вывода

Если ваша функция возвращает массив Array или какой-либо другой сложный тип, ей может потребоваться выделить память. К сожалению, выделение памяти и обратное действие — сборка мусора — часто являются существенными узкими местами.

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

julia> function xinc(x)
           return [x, x+1, x+2]
       end;

julia> function loopinc()
           y = 0
           for i = 1:10^7
               ret = xinc(i)
               y += ret[2]
           end
           return y
       end;

julia> function xinc!(ret::AbstractVector{T}, x::T) where T
           ret[1] = x
           ret[2] = x+1
           ret[3] = x+2
           nothing
       end;

julia> function loopinc_prealloc()
           ret = Vector{Int}(undef, 3)
           y = 0
           for i = 1:10^7
               xinc!(ret, i)
               y += ret[2]
           end
           return y
       end;

Результаты измерений времени:

julia> @time loopinc()
  0.529894 seconds (40.00 M allocations: 1.490 GiB, 12.14% gc time)
50000015000000

julia> @time loopinc_prealloc()
  0.030850 seconds (6 allocations: 288 bytes)
50000015000000

END_OF_DOCUMENT_MARKER

Предварительное выделение памяти имеет другие преимущества, например, позволяя вызывающей стороне управлять типом «вывода» алгоритма. В приведенном выше примере мы могли бы передать SubArray вместо Array, если бы захотели.

Доведенное до крайности, предварительное выделение памяти может ухудшить читаемость вашего кода, поэтому могут потребоваться измерения производительности и некоторый здравый смысл. Однако для «векторизованных» (элементных) функций может быть удобен синтаксис x .= f.(y) для операций «на месте» с объединенными циклами и без временных массивов (см. синтаксис точки для векторизации функций).

Используйте `MutableArithmetics` для большего контроля над выделением для изменяемых типов арифметики

Некоторые подтипы Number, такие как BigInt или BigFloat, могут быть реализованы как типы mutable struct, или они могут иметь изменяемые компоненты. Арифметические интерфейсы в Julia Base обычно отдают предпочтение удобству перед эффективностью в таких случаях, поэтому их использование простым способом может привести к неэффективной производительности. С другой стороны, абстракции пакета MutableArithmetics позволяют использовать изменчивость таких типов для написания быстрого кода, выделяющего только необходимый объем памяти. MutableArithmetics также позволяет явно копировать значения изменяемых типов арифметики при необходимости. MutableArithmetics - это пользовательский пакет и не связан с проектом Julia.

Больше точек: Объединить векторизованные операции

В Julia есть специальный синтаксис точки, который преобразует любую скалярную функцию в вызов «векторизованной» функции, а любой оператор — в «векторизованный» оператор со специальным свойством, что вложенные «вызовы с точкой» объединяются: они объединяются на уровне синтаксиса в один цикл без выделения временных массивов. Если вы используете .= и подобные операторы присваивания, результат также можно сохранить на месте в предварительно выделенном массиве (см. выше).

В контексте линейной алгебры это означает, что даже если операции, такие как vector + vector и vector * scalar, определены, может быть выгодно использовать vector .+ vector и vector .* scalar, так как полученные циклы можно объединить с окружающими вычислениями. Например, рассмотрим две функции:

julia> f(x) = 3x.^2 + 4x + 7x.^3;

julia> fdot(x) = @. 3x^2 + 4x + 7x^3; # equivalent to 3 .* x.^2 .+ 4 .* x .+ 7 .* x.^3

И f, и fdot вычисляют одно и то же. Однако, fdot (определённая с помощью макроса @.) значительно быстрее при применении к массиву:

julia> x = rand(10^6);

julia> @time f(x);
  0.019049 seconds (16 allocations: 45.777 MiB, 18.59% gc time)

julia> @time fdot(x);
  0.002790 seconds (6 allocations: 7.630 MiB)

julia> @time f.(x);
  0.002626 seconds (8 allocations: 7.630 MiB)

То есть, fdot(x) в десять раз быстрее и выделяет в 1/6 памяти, чем f(x), потому что каждая операция * и + в f(x) выделяет новый временный массив и выполняется в отдельном цикле. В этом примере f.(x) так же быстра, как fdot(x), но во многих контекстах более удобно разбросать несколько точек в ваших выражениях, чем определять отдельную функцию для каждой векторизованной операции.

Меньше точек: Разделить определенные промежуточные широковещательные передачи

Объединение циклов с точкой, упомянутое выше, позволяет кратко и идиоматично выразить высокопроизводительные операции. Однако важно помнить, что объединённая операция будет вычисляться на каждой итерации широковещательной передачи. Это означает, что в некоторых ситуациях, особенно при наличии составных или многомерных широковещательных передач, выражение с вызовами точек может вычислять функцию чаще, чем ожидается. Предположим, что мы хотим создать случайную матрицу, у которой строки имеют евклидову норму один. Мы можем написать что-то вроде следующего:

julia> x = rand(1000, 1000);

julia> d = sum(abs2, x; dims=2);

julia> @time x ./= sqrt.(d);
  0.002049 seconds (4 allocations: 96 bytes)

Это сработает. Однако это выражение фактически повторно вычислит sqrt(d[i]) для каждого элемента в строке x[i, :], что означает, что вычисляется гораздо больше квадратных корней, чем необходимо. Чтобы увидеть, по каким именно индексам будет происходить широковещательная передача, мы можем вызвать Broadcast.combine_axes для аргументов объединённого выражения. Это вернёт кортеж диапазонов, записи которого соответствуют осям итерации; произведение длин этих диапазонов будет общим числом вызовов объединённой операции.

Следовательно, когда некоторые компоненты выражения широковещательной передачи являются постоянными вдоль оси, как sqrt вдоль второй оси в предыдущем примере, существует возможность улучшить производительность, принудительно «разделяя» эти компоненты, т.е. предварительно выделяя результат широковещательной операции и повторно используя кэшированное значение вдоль его постоянной оси. Некоторые такие потенциальные подходы — использовать временные переменные, обернуть компоненты выражения с точкой в identity, или использовать эквивалентную собственно векторизованную (но не объединённую) функцию.

julia> @time let s = sqrt.(d); x ./= s end;
  0.000809 seconds (5 allocations: 8.031 KiB)

julia> @time x ./= identity(sqrt.(d));
  0.000608 seconds (5 allocations: 8.031 KiB)

julia> @time x ./= map(sqrt, d);
  0.000611 seconds (4 allocations: 8.016 KiB)

Любой из этих вариантов даёт примерно трехкратное ускорение за счёт выделения памяти; для больших широковещательных данных это ускорение может быть асимптотически очень большим.

Рассмотрите использование представлений для срезов

В Julia выражение «срез» массива, такое как array[1:5, :] создаёт копию этих данных (за исключением левой части присваивания, где array[1:5, :] = ... выполняет присваивание на месте в эту часть array). Если вы выполняете много операций над срезoм, это может быть выгодно для производительности, так как работать с меньшей непрерывной копией эффективнее, чем индексировать исходный массив. С другой стороны, если вы выполняете только несколько простых операций над срезoм, затраты на выделение памяти и копирование могут быть существенными.

Альтернативой является создание «представления» массива, которое является объектом массива (SubArray) , фактически ссылающимся на данные исходного массива на месте без копирования. (Если вы записываете в представление, оно также изменяет данные исходного массива.) Это можно сделать для отдельных срезов, вызвав view, или проще для целого выражения или блока кода, поместив @views перед этим выражением. Например:

julia> fcopy(x) = sum(x[2:end-1]);

julia> @views fview(x) = sum(x[2:end-1]);

julia> x = rand(10^6);

julia> @time fcopy(x);
  0.003051 seconds (3 allocations: 7.629 MB)

julia> @time fview(x);
  0.001020 seconds (1 allocation: 16 bytes)

Обратите внимание на трехкратное ускорение и уменьшение выделения памяти версии функции fview.

Копирование данных не всегда плохо

Массивы хранятся непрерывно в памяти, что способствует векторизации CPU и меньшему числу обращений к памяти за счёт кэширования. По тем же причинам рекомендуется обращаться к массивам в порядке следования столбцов (см. выше). Неправильные шаблоны доступа и несмежные представления могут значительно замедлить вычисления с массивами из-за не последовательных обращений к памяти.

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

julia> using Random

julia> A = randn(3000, 3000);

julia> x = randn(2000);

julia> inds = shuffle(1:3000)[1:2000];

julia> function iterated_neural_network(A, x, depth)
           for _ in 1:depth
               x .= max.(0, A * x)
           end
           argmax(x)
       end

julia> @time iterated_neural_network(view(A, inds, inds), x, 10)
  0.324903 seconds (12 allocations: 157.562 KiB)
1569

julia> @time iterated_neural_network(A[inds, inds], x, 10)
  0.054576 seconds (13 allocations: 30.671 MiB, 13.33% gc time)
1569

При наличии достаточного объёма памяти стоимость копирования представления в массив перевешивается улучшением скорости многократного умножения матриц на непрерывном массиве.

Рассмотрите StaticArrays.jl для небольших операций с фиксированными векторами/матрицами

Если ваше приложение включает в себя множество небольших (< 100 элементов) массивов фиксированного размера (т. е. размер известен до выполнения), то вы можете рассмотреть использование пакета StaticArrays.jl. Этот пакет позволяет представлять такие массивы таким образом, чтобы избежать ненужных выделений памяти и позволяет компилятору специализировать код для размера массива, например, путём полного развёртывания векторных операций (удаление циклов) и хранения элементов в регистрах процессора.

Например, если вы выполняете вычисления с 2D геометриями, вы можете выполнять много вычислений с векторами с 2 компонентами. Используя тип SVector из StaticArrays.jl, вы можете использовать удобную векторную нотацию и операции, такие как norm(3v - w) над векторами v и w, позволяя компилятору свернуть код до минимального вычисления, эквивалентного @inbounds hypot(3v[1]-w[1], 3v[2]-w[2]).

Избегайте интерполяции строк для ввода/вывода

При записи данных в файл (или другой устройстве ввода/вывода) формирование дополнительных промежуточных строк — источник накладных расходов. Вместо:

println(file, "$a $b")

используйте:

println(file, a, " ", b)

В первой версии кода создаётся строка, а затем она записывается в файл, а во второй версии значения записываются непосредственно в файл. Также обратите внимание, что в некоторых случаях интерполяция строк может быть сложнее для чтения. Рассмотрите:

println(file, "$(f(a))$(f(b))")

по сравнению с:

println(file, f(a), f(b))

Оптимизируйте сетевой ввод/вывод во время параллельного выполнения

При выполнении удалённой функции параллельно:

using Distributed

responses = Vector{Any}(undef, nworkers())
@sync begin
    for (idx, pid) in enumerate(workers())
        @async responses[idx] = remotecall_fetch(foo, pid, args...)
    end
end

быстрее, чем:

using Distributed

refs = Vector{Any}(undef, nworkers())
for (idx, pid) in enumerate(workers())
    refs[idx] = @spawnat pid foo(args...)
end
responses = [fetch(r) for r in refs]

Первый вариант приводит к одному сетевому обмену с каждым рабочим узлом, в то время как второй вариант приводит к двум сетевым вызовам — сначала от @spawnat, а затем из-за fetch (или даже wait). fetch/wait также выполняется последовательно, что приводит к более низкой производительности в целом.

Исправьте предупреждения о устаревании

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

Настройки

Вот некоторые незначительные моменты, которые могут помочь в узких внутренних циклах.

Избегайте ненужных массивов. Например, вместо sum([x,y,z]) используйте x+y+z.
Используйте abs2(z) вместо abs(z)^2 для сложных z. В общем, старайтесь переписывать код, чтобы использовать abs2 вместо abs для сложных аргументов.
Используйте div(x,y) для усечения целочисленного деления вместо trunc(x/y), fld(x,y) вместо floor(x/y) и cld(x,y) вместо ceil(x/y).

Аннотации производительности

Иногда вы можете включить лучшую оптимизацию, гарантировав определенные свойства программы.

Используйте @inbounds, чтобы исключить проверку границ массива в выражениях. Убедитесь в этом. Если индексы когда-либо выходят за пределы границ, вы можете столкнуться с ошибками или скрытым повреждением.
Используйте @fastmath, чтобы разрешить оптимизации с плавающей запятой, которые верны для вещественных чисел, но приводят к различиям для чисел IEEE. Будьте осторожны при этом, так как это может изменить численные результаты. Это соответствует опции -ffast-math clang.
Напишите @simd перед for циклами, чтобы гарантировать, что итерации независимы и могут быть переупорядочены. Обратите внимание, что во многих случаях Julia может автоматически векторизовать код без макроса @simd; это полезно только в тех случаях, когда такая трансформация в противном случае была бы незаконной, включая случаи, когда разрешается повторная ассоциативность чисел с плавающей запятой и игнорируются зависимые обращения к памяти (@simd ivdep). Снова будьте очень осторожны при утверждении @simd, так как неверная аннотация цикла с зависимыми итерациями может привести к неожиданным результатам. В частности, обратите внимание, что setindex! для некоторых AbstractArray подтипов изначально зависит от порядка итераций. Эта функция экспериментальна и может измениться или исчезнуть в будущих версиях Julia.

Общий прием использования 1:n для индексации в AbstractArray не безопасен, если массив использует нестандартную индексацию, и может вызвать ошибку сегментации, если проверка границ отключена. Используйте LinearIndices(x) или eachindex(x) вместо этого (см. также Массивы с пользовательскими индексами).

Хотя @simd необходимо размещать непосредственно перед вложенным for циклом, @inbounds и @fastmath могут применяться как к отдельным выражениям, так и ко всем выражениям, которые появляются внутри вложенных блоков кода, например, используя @inbounds begin или @inbounds for ....

Вот пример с аннотациями @inbounds и @simd (здесь мы используем @noinline для предотвращения попыток оптимизатора быть слишком умным и нарушать наш бенчмарк):

@noinline function inner(x, y)
    s = zero(eltype(x))
    for i=eachindex(x)
        @inbounds s += x[i]*y[i]
    end
    return s
end

@noinline function innersimd(x, y)
    s = zero(eltype(x))
    @simd for i = eachindex(x)
        @inbounds s += x[i] * y[i]
    end
    return s
end

function timeit(n, reps)
    x = rand(Float32, n)
    y = rand(Float32, n)
    s = zero(Float64)
    time = @elapsed for j in 1:reps
        s += inner(x, y)
    end
    println("GFlop/sec        = ", 2n*reps / time*1E-9)
    time = @elapsed for j in 1:reps
        s += innersimd(x, y)
    end
    println("GFlop/sec (SIMD) = ", 2n*reps / time*1E-9)
end

timeit(1000, 1000)

На компьютере с процессором Intel Core i5 2,4 ГГц это даёт:

GFlop/sec        = 1.9467069505224963
GFlop/sec (SIMD) = 17.578554163920018

(GFlop/sec измеряет производительность, и большие числа лучше.)

Вот пример со всеми тремя видами аннотаций. Эта программа сначала вычисляет конечную разность одномерного массива, а затем оценивает L2-норму результата:

function init!(u::Vector)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds @simd for i in 1:n #by asserting that `u` is a `Vector` we can assume it has 1-based indexing
        u[i] = sin(2pi*dx*i)
    end
end

function deriv!(u::Vector, du)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds du[1] = (u[2] - u[1]) / dx
    @fastmath @inbounds @simd for i in 2:n-1
        du[i] = (u[i+1] - u[i-1]) / (2*dx)
    end
    @fastmath @inbounds du[n] = (u[n] - u[n-1]) / dx
end

function mynorm(u::Vector)
    n = length(u)
    T = eltype(u)
    s = zero(T)
    @fastmath @inbounds @simd for i in 1:n
        s += u[i]^2
    end
    @fastmath @inbounds return sqrt(s)
end

function main()
    n = 2000
    u = Vector{Float64}(undef, n)
    init!(u)
    du = similar(u)

    deriv!(u, du)
    nu = mynorm(du)

    @time for i in 1:10^6
        deriv!(u, du)
        nu = mynorm(du)
    end

    println(nu)
end

main()

На компьютере с процессором Intel Core i7 2,7 ГГц это даёт:

$ julia wave.jl;
  1.207814709 seconds
4.443986180758249

$ julia --math-mode=ieee wave.jl;
  4.487083643 seconds
4.443986180758249

Здесь опция --math-mode=ieee отключает макрос @fastmath, чтобы мы могли сравнить результаты.

В этом случае ускорение благодаря @fastmath составляет примерно 3,7 раза. Это необычно большое значение — обычно ускорение будет меньше. (В этом конкретном примере рабочая область бенчмарка достаточно мала, чтобы поместиться в кэш L1 процессора, поэтому задержка доступа к памяти не играет роли, а время вычислений определяется использованием ЦП. В многих реальных программах это не так.) Кроме того, в этом случае эта оптимизация не изменяет результат — обычно результат будет немного отличаться. В некоторых случаях, особенно для численных неустойчивых алгоритмов, результат может сильно отличаться.

Аннотация @fastmath переупорядочивает выражения с плавающей запятой, например, изменяя порядок вычислений или предполагая, что некоторые особые случаи (inf, nan) не могут произойти. В этом случае (и на этом конкретном компьютере) основное различие заключается в том, что выражение 1 / (2*dx) в функции deriv выносится за цикл (т.е. вычисляется вне цикла), как если бы было написано idx = 1 / (2*dx). В цикле выражение ... / (2*dx) тогда становится ... * idx, что намного быстрее для вычисления. Конечно, как сама оптимизация, применяемая компилятором, так и полученное ускорение очень зависят от аппаратного обеспечения. Вы можете изучить изменения в сгенерированном коде, используя функцию Julia code_native.

Обратите внимание, что @fastmath также предполагает, что NaN не будут происходить во время вычислений, что может привести к неожиданному поведению:

julia> f(x) = isnan(x);

julia> f(NaN)
true

julia> f_fast(x) = @fastmath isnan(x);

julia> f_fast(NaN)
false

Обработка чисел с пониженной точностью как нулей

Числа с пониженной точностью, ранее называвшиеся числами с пониженной точностью, полезны во многих контекстах, но вносят штраф за производительность на некоторых аппаратных платформах. Вызов set_zero_subnormals(true) даёт разрешение операциям с плавающей запятой обрабатывать поднормальные входные или выходные данные как нули, что может улучшить производительность на некоторых аппаратных платформах. Вызов set_zero_subnormals(false) обеспечивает строгое поведение IEEE для чисел с пониженной точностью.

Ниже приведён пример, где числа с пониженной точностью заметно влияют на производительность на некоторых аппаратных платформах:

function timestep(b::Vector{T}, a::Vector{T}, Δt::T) where T
    @assert length(a)==length(b)
    n = length(b)
    b[1] = 1                            # Boundary condition
    for i=2:n-1
        b[i] = a[i] + (a[i-1] - T(2)*a[i] + a[i+1]) * Δt
    end
    b[n] = 0                            # Boundary condition
end

function heatflow(a::Vector{T}, nstep::Integer) where T
    b = similar(a)
    for t=1:div(nstep,2)                # Assume nstep is even
        timestep(b,a,T(0.1))
        timestep(a,b,T(0.1))
    end
end

heatflow(zeros(Float32,10),2)           # Force compilation
for trial=1:6
    a = zeros(Float32,1000)
    set_zero_subnormals(iseven(trial))  # Odd trials use strict IEEE arithmetic
    @time heatflow(a,1000)
end

Это даёт вывод, похожий на

  0.002202 seconds (1 allocation: 4.063 KiB)
  0.001502 seconds (1 allocation: 4.063 KiB)
  0.002139 seconds (1 allocation: 4.063 KiB)
  0.001454 seconds (1 allocation: 4.063 KiB)
  0.002115 seconds (1 allocation: 4.063 KiB)
  0.001455 seconds (1 allocation: 4.063 KiB)

Обратите внимание, как каждая чётная итерация значительно быстрее.

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

Обработка чисел с пониженной точностью как нулей должна использоваться с осторожностью, поскольку это нарушает некоторые тождества, например, x-y == 0 подразумевает x == y.

julia> x = 3f-38; y = 2f-38;

julia> set_zero_subnormals(true); (x - y, x == y)
(0.0f0, false)

julia> set_zero_subnormals(false); (x - y, x == y)
(1.0000001f-38, false)

В некоторых приложениях альтернативой нулевым числам с пониженной точностью является введение небольшого шума. Например, вместо инициализации a нулями, инициализируйте её:

a = rand(Float32,1000) * 1.f-9

`@code_warntype`

Макрос @code_warntype (или его функциональный вариант code_warntype) иногда бывает полезен для диагностики проблем, связанных с типами. Вот пример:

julia> @noinline pos(x) = x < 0 ? 0 : x;

julia> function f(x)
           y = pos(x)
           return sin(y*x + 1)
       end;

julia> @code_warntype f(3.2)
MethodInstance for f(::Float64)
  from f(x) @ Main REPL[9]:1
Arguments
  #self#::Core.Const(f)
  x::Float64
Locals
  y::Union{Float64, Int64}
Body::Float64
1 ─      (y = Main.pos(x))
│   %2 = (y * x)::Float64
│   %3 = (%2 + 1)::Float64
│   %4 = Main.sin(%3)::Float64
└──      return %4

Интерпретация вывода @code_warntype, как и его аналогов @code_lowered, @code_typed, @code_llvm и @code_native, требует некоторой практики. Ваш код представлен в виде, который был сильно переработан на пути к генерации скомпилированного машинного кода. Большинство выражений аннотированы типом, указанным в ::T (например, T может быть Float64). Наиболее важной характеристикой @code_warntype является то, что неконкретные типы отображаются красным цветом; поскольку этот документ написан на языке разметки Markdown, в котором нет цвета, в этом документе красные тексты обозначены заглавными буквами.

Вверху показан выведенный возвращаемый тип функции как Body::Float64. Следующие строки представляют тело f в форме SSA-IR Julia. Номерные блоки — метки и представляют цели для переходов (через goto) в вашем коде. Рассматривая тело, вы можете увидеть, что в первую очередь вызывается pos, а возвращаемое значение определено как тип Union Union{Float64, Int64}, отображаемый заглавными буквами, поскольку он является неконкретным типом. Это означает, что мы не можем точно определить возвращаемый тип pos на основе входных типов. Однако результат y*x — это Float64, независимо от того, является ли y Float64 или Int64. Итоговым результатом является то, что f(x::Float64) не будет иметь нестабильный тип вывода, даже если некоторые промежуточные вычисления имеют нестабильный тип.

Как вы будете использовать эту информацию — решать вам. Очевидно, было бы намного лучше исправить pos так, чтобы он был устойчив к типу: если вы это сделаете, все переменные в f будут конкретными, и его производительность будет оптимальной. Однако есть обстоятельства, когда такая временная неустойчивость типа может не иметь большого значения: например, если pos никогда не используется изолированно, тот факт, что выход f устойчив к типу (для входных данных Float64) защитит последующий код от распространяющегося влияния неустойчивости типа. Это особенно актуально в тех случаях, когда исправление неустойчивости типа затруднено или невозможно. В таких случаях приведенные выше советы (например, добавление аннотаций типа и/или разделение функций) являются лучшими инструментами для ограничения «повреждений», вызванных неустойчивостью типа. Также обратите внимание, что даже в Julia Base есть функции, которые неустойчивы к типу. Например, функция findfirst возвращает индекс в массиве, где найден ключ, или nothing если он не найден, что является явной неустойчивостью типа. Чтобы легче было найти неустойчивости типа, которые, вероятно, важны, Union , содержащие либо missing , либо nothing , выделяются жёлтым цветом вместо красного.

Следующие примеры могут помочь вам интерпретировать выражения, помеченные как содержащие типы, не являющиеся листьями:

Тело функции, начинающееся с Body::Union{T1,T2})
- Интерпретация: функция с неустойчивым типом возвращаемого значения
- Предложение: сделайте тип возвращаемого значения устойчивым к типу, даже если вам придётся его аннотировать
invoke Main.g(%%x::Int64)::Union{Float64, Int64}
- Интерпретация: вызов функции с неустойчивым типом g.
- Предложение: исправить функцию или, если необходимо, аннотировать возвращаемое значение
invoke Base.getindex(%%x::Array{Any,1}, 1::Int64)::Any
- Интерпретация: доступ к элементам массивов с нечётко определённым типом
- Предложение: используйте массивы с лучше определёнными типами или, при необходимости, аннотируйте тип отдельных обращений к элементам
Base.getfield(%%x, :(:data))::Array{Float64,N} where N
- Интерпретация: получение поля, которое имеет тип, не являющийся листом. В данном случае тип x, скажем, ArrayContainer, имел поле data::Array{T}. Но Array также требует размерность N, чтобы быть конкретным типом.
- Предложение: используйте конкретные типы, такие как Array{T,3} или Array{T,N}, где N теперь является параметром ArrayContainer

Производительность захваченной переменной

Рассмотрим следующий пример, который определяет внутреннюю функцию:

function abmult(r::Int)
    if r < 0
        r = -r
    end
    f = x -> x * r
    return f
end

Функция abmult возвращает функцию f , которая умножает свой аргумент на абсолютное значение r . Внутренняя функция, назначенная f , называется «замыканием». Внутренние функции также используются языком для do-блоков и для выражений генераторов.

Этот стиль кода представляет собой проблему производительности для языка. Парсер, переводя его в инструкции более низкого уровня, существенно перестраивает приведенный выше код, вынося внутреннюю функцию в отдельный блок кода. «Захваченные» переменные, такие как r , которые используются внутренними функциями и их охватывающей областью видимости, также переносятся в выделенный на куче «ящик» (box), доступный как внутренним, так и внешним функциям, поскольку язык определяет, что r во внутренней области видимости должно быть идентично r во внешней области видимости даже после того, как внешняя область видимости (или другая внутренняя функция) изменит r.

Обсуждение в предыдущем абзаце касалось «парсера», то есть этапа компиляции, который выполняется при первом загрузке модуля, содержащего abmult , в отличие от последующего этапа, когда он впервые вызывается. Парсер не «знает», что Int — это фиксированный тип или что оператор r = -r преобразует Int в другой Int . Магия вывода типов происходит на более позднем этапе компиляции.

Таким образом, парсер не знает, что r имеет фиксированный тип (Int), ни то, что r не изменяет значение после создания внутренней функции (что делает ящик (box) ненужным). Поэтому парсер генерирует код для ящика, который содержит объект с абстрактным типом, таким как Any , что требует диспетчеризации типа во время выполнения для каждого появления r . Это можно проверить, применив @code_warntype к приведенной выше функции. Как создание ящика, так и диспетчеризация типа во время выполнения могут привести к потере производительности.

Если захваченные переменные используются в критичной для производительности части кода, следующие советы помогают обеспечить, чтобы их использование было эффективным. Во-первых, если известно, что захваченная переменная не изменяет свой тип, это можно явно указать с помощью аннотации типа (к переменной, а не к правой части):

function abmult2(r0::Int)
    r::Int = r0
    if r < 0
        r = -r
    end
    f = x -> x * r
    return f
end

Аннотация типа частично восстанавливает производительность, потерянную из-за захвата, поскольку парсер может связать конкретный тип с объектом в ящике. Далее, если захваченной переменной вообще не нужно создавать ящик (потому что она не будет переприсваиваться после создания замыкания), это можно указать с помощью let блоков следующим образом.

function abmult3(r::Int)
    if r < 0
        r = -r
    end
    f = let r = r
            x -> x * r
    end
    return f
end

Блок let создаёт новую переменную r , область видимости которой ограничена только внутренней функцией. Второй приём восстанавливает полную производительность языка при наличии захваченных переменных. Обратите внимание, что это быстро развивающийся аспект компилятора, и в будущих выпусках, вероятно, не потребуется такая степень аннотации программистом для достижения производительности. Тем временем некоторые пользовательские пакеты, такие как FastClosures, автоматизируют вставку let операторов, как в abmult3.

Многопоточность и линейная алгебра

Этот раздел относится к многопоточному коду Julia, который в каждом потоке выполняет операции линейной алгебры. Действительно, эти операции линейной алгебры включают вызовы BLAS/LAPACK, которые сами по себе многопоточны. В этом случае необходимо убедиться, что ядра не перегружены из-за двух разных типов многопоточности.

Julia компилирует и использует свою собственную копию OpenBLAS для линейной алгебры, количество потоков которой контролируется переменной среды OPENBLAS_NUM_THREADS. Её можно установить как параметр командной строки при запуске Julia или изменить во время сессии Julia с помощью BLAS.set_num_threads(N) (подмодуль BLAS экспортируется using LinearAlgebra). Его текущее значение можно получить с помощью BLAS.get_num_threads().

Если пользователь ничего не указывает, Julia пытается выбрать разумное значение для количества потоков OpenBLAS (например, на основе платформы, версии Julia и т. д.). Однако, как правило, рекомендуется проверять и устанавливать значение вручную. Поведение OpenBLAS таково:

Если OPENBLAS_NUM_THREADS=1, OpenBLAS использует вызывающий(ие) поток(и) Julia, то есть «живёт» в потоке Julia, который выполняет вычисления.
Если OPENBLAS_NUM_THREADS=N>1, OpenBLAS создаёт и управляет своим собственным пулом потоков (N в целом). Существует только один пул потоков OpenBLAS, используемый всеми потоками Julia.

При запуске Julia в многопоточном режиме с JULIA_NUM_THREADS=X рекомендуется установить OPENBLAS_NUM_THREADS=1 . Учитывая описанное выше поведение, увеличение количества потоков BLAS до N>1 очень легко может привести к ухудшению производительности, особенно когда N<<X. Однако это всего лишь эмпирическое правило, и наилучший способ установить каждое количество потоков — это эксперимент на вашей конкретной программе.

Альтернативные бэкэнды линейной алгебры

В качестве альтернативы OpenBLAS существует несколько других бэкэндов, которые могут помочь в повышении производительности линейной алгебры. Яркими примерами являются MKL.jl и AppleAccelerate.jl.

Это внешние пакеты, поэтому мы не будем подробно их обсуждать здесь. Обратитесь к соответствующей документации (особенно потому, что они ведут себя по-разному, чем OpenBLAS, в отношении многопоточности).

Задержка выполнения, загрузка пакетов и время предварительной компиляции пакетов

Сокращение времени до первого графика и т. д.

В первый раз при вызове метода julia (и всех вызываемых им методах или методах, которые могут быть статически определены) он будет компилироваться. Семья макросов @time иллюстрирует это.

julia> foo() = rand(2,2) * rand(2,2)
foo (generic function with 1 method)

julia> @time @eval foo();
  0.252395 seconds (1.12 M allocations: 56.178 MiB, 2.93% gc time, 98.12% compilation time)

julia> @time @eval foo();
  0.000156 seconds (63 allocations: 2.453 KiB)

Обратите внимание, что @time @eval лучше для измерения времени компиляции, потому что без @eval часть компиляции может быть уже выполнена до начала измерения времени.

При разработке пакета вы можете улучшить опыт своих пользователей с помощью предварительной компиляции, чтобы при использовании пакета код, который они используют, был уже скомпилирован. Для эффективной предварительной компиляции кода пакета рекомендуется использовать PrecompileTools.jl для выполнения «задачи предварительной компиляции» во время предварительной компиляции, которая является представительной для типичного использования пакета, что позволит кешировать нативный скомпилированный код в кэш пакета pkgimage , значительно сокращая время до первого выполнения (часто обозначаемое как TTFX) для такого использования.

Обратите внимание, что задачи предварительной компиляции PrecompileTools.jl можно отключить и иногда настроить через «Настройки», если вы не хотите тратить дополнительное время на предварительную компиляцию, что может быть актуально во время разработки пакета.

Сокращение времени загрузки пакета

Уменьшение времени, необходимого для загрузки пакета, обычно полезно. Общие рекомендации для разработчиков пакетов включают:

Сократите зависимости до необходимых. Рассмотрите использование расширений пакетов, чтобы обеспечить совместимость с другими пакетами, не увеличивая объем основных зависимостей.
Избегайте использования __init__() функций, если есть альтернатива, особенно тех, которые могут вызывать много компиляций или долго выполняться.
По возможности исправьте несоответствия среди зависимостей и в коде вашего пакета.

Инструмент @time_imports может быть полезен в REPL для анализа вышеперечисленных факторов.

julia> @time @time_imports using Plots
      0.5 ms  Printf
     16.4 ms  Dates
      0.7 ms  Statistics
               ┌ 23.8 ms SuiteSparse_jll.__init__() 86.11% compilation time (100% recompilation)
     90.1 ms  SuiteSparse_jll 91.57% compilation time (82% recompilation)
      0.9 ms  Serialization
               ┌ 39.8 ms SparseArrays.CHOLMOD.__init__() 99.47% compilation time (100% recompilation)
    166.9 ms  SparseArrays 23.74% compilation time (100% recompilation)
      0.4 ms  Statistics → SparseArraysExt
      0.5 ms  TOML
      8.0 ms  Preferences
      0.3 ms  PrecompileTools
      0.2 ms  Reexport
... many deps omitted for example ...
      1.4 ms  Tar
               ┌ 73.8 ms p7zip_jll.__init__() 99.93% compilation time (100% recompilation)
     79.4 ms  p7zip_jll 92.91% compilation time (100% recompilation)
               ┌ 27.7 ms GR.GRPreferences.__init__() 99.77% compilation time (100% recompilation)
     43.0 ms  GR 64.26% compilation time (100% recompilation)
               ┌ 2.1 ms Plots.__init__() 91.80% compilation time (100% recompilation)
    300.9 ms  Plots 0.65% compilation time (100% recompilation)
  1.795602 seconds (3.33 M allocations: 190.153 MiB, 7.91% gc time, 39.45% compilation time: 97% of which was recompilation)

Обратите внимание, что в этом примере загружено несколько пакетов, некоторые с __init__() функциями, некоторые из которых вызывают компиляцию, включая повторную компиляцию. Повторная компиляция вызывается, когда предыдущие пакеты аннулируют методы, а затем, когда следующие пакеты выполняют свою __init__() функцию, происходит повторная компиляция, прежде чем код может быть выполнен.

Кроме того, обратите внимание на Statistics расширение SparseArraysExt, которое активировано, потому что SparseArrays присутствует в дереве зависимостей. То есть, см. 0.4 ms Statistics → SparseArraysExt.

Этот отчет предоставляет хорошую возможность оценить, оправдывает ли стоимость времени загрузки зависимостей предоставляемую функциональность. Также Pkg утилита why может быть использована для вывода причины существования косвенной зависимости.

(CustomPackage) pkg> why FFMPEG_jll
  Plots → FFMPEG → FFMPEG_jll
  Plots → GR → GR_jll → FFMPEG_jll

Или чтобы увидеть косвенные зависимости, которые вносит пакет, вы можете pkg> rm пакет, увидеть зависимости, которые удалены из манифеста, а затем восстановить изменения с помощью pkg> undo.

Если время загрузки доминируется медленными __init__() методами с компиляцией, то один из способов подробной идентификации компилируемого кода заключается в использовании аргументов Julia --trace-compile=stderr, которые будут сообщать о precompile операторе каждый раз, когда компилируется метод. Например, полная настройка будет:

$ julia --startup-file=no --trace-compile=stderr
julia> @time @time_imports using CustomPackage
...

Обратите внимание на --startup-file=no, который помогает изолировать тест от пакетов, которые могут быть в вашем startup.jl.

Более подробный анализ причин повторной компиляции можно получить с помощью пакета SnoopCompile.

© 2009–2024 Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors
Licensed under the MIT License.
https://docs.julialang.org/en/v1.11/manual/performance-tips/

Рекомендации по производительности

Измеряйте производительность с помощью @time и обращайте внимание на выделение памяти

Измеряйте производительность с помощью `@time` и обращайте внимание на выделение памяти