Советы по производительности

В следующих разделах мы кратко рассмотрим несколько техник, которые помогут сделать ваш код на 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, строк или «нестабильных по типу» объектов, тип которых известен только во время выполнения). Выделение (или освобождение) таких блоков памяти может потребовать дорогостоящего системного вызова (например, через malloc в C), и их необходимо отслеживать для сборки мусора. В отличие от неизменяемых значений, таких как числа (за исключением больших чисел), кортежи и неизменяемые structs, которые можно хранить намного дешевле, например, в стеке или регистре ЦП, поэтому, как правило, не нужно беспокоиться о производительности «выделения».

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

В этом конкретном случае выделение памяти связано с использованием нестабильной по типу глобальной переменной 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 поможет вам визуализировать результаты профилирования.
Пакет Traceur может помочь вам найти распространенные проблемы с производительностью в вашем коде.
Неожиданно большие выделения памяти — как сообщается @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.

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

Такие же лучшие практики также работают для типов контейнеров:

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 учитывает структуру памяти массива Matrix по столбцам и заполняет его по одному столбцу за раз. Кроме того, 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

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

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

Ещё точки: Объединение векторизованных операций

У 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) в десять раз быстрее и выделяет в 6 раз меньше памяти, чем f(x), потому что каждая операция * и + в f(x) выделяет новый временный массив и выполняется в отдельном цикле. В этом примере f.(x) так же быстра, как fdot(x), но во многих контекстах удобнее расставлять точки в ваших выражениях, чем определять отдельную функцию для каждой векторизованной операции.

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

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

Альтернативой является создание «представления» массива, которое является объектом массива (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)

Обратите внимание на 3-кратное увеличение скорости и уменьшение выделения памяти версии 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 , которые используются как внутренней, так и окружающей функциями, также извлекаются в «ящик» в куче, доступный как внутренним, так и внешним функциям, потому что язык требует, чтобы r во внутреннем пространстве имён было идентично r во внешнем пространстве имён даже после того, как внешнее пространство имён (или другая внутренняя функция) изменяет r.

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

Таким образом, парсер не знает, что r имеет фиксированный тип (Int), или что r не изменяет значение после создания внутренней функции (так что «ящик» не нужен). Поэтому парсер генерирует код для ящика, который содержит объект с абстрактным типом, таким как 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 в отношении многопоточности).

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

Советы по производительности

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

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