Performance Tips

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

Performance critical code should be inside a function

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

Использование функций важно не только для производительности: функции более переиспользуемы и тестируемы, а также проясняют, какие шаги выполняются и каковы их входные и выходные данные, Write functions, not just scripts также является рекомендацией Стиля Julia.

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

Avoid untyped global variables

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

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

const DEFAULT_VAL = 0

Если глобальная переменная известна как всегда имеющая один и тот же тип, the type should be annotated.

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

global x = rand(1000)

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

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

Note

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

В следующей сессии REPL:

julia> x = 1.0

равно:

julia> global x = 1.0

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

Measure performance with `@time` and pay attention to memory allocation

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

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

Note

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

Tools

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

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

Avoid containers with abstract type parameters

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

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

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, их можно будет хранить более эффективно:

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}.

См. также обсуждение под Parametric Types.

Type declarations

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

Avoid fields with abstract type

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

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.

Avoid fields with abstract containers

Те же лучшие практики также применимы к типам контейнеров:

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

Annotate values taken from untyped locations

Часто удобно работать с структурами данных, которые могут содержать значения любого типа (массивы типа 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 ухудшает производительность. Чтобы написать производительный код, включающий типы, создаваемые во время выполнения, используйте function-barrier technique, обсужденный ниже, и убедитесь, что созданный тип присутствует среди типов аргументов функции ядра, чтобы операции ядра были правильно специализированы компилятором. Например, в приведенном выше фрагменте, как только b создан, его можно передать другой функции k, ядру. Если, например, функция k объявляет b как аргумент типа Complex{T}, где T — это параметр типа, тогда аннотация типа, appearing in an assignment statement within k of the form:

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

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

Be aware of when Julia avoids specializing

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

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

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 обычно не специализировала бы этот вызов метода. Вам нужно проверить method internals, если вы хотите увидеть, генерируются ли специализации, когда изменяются типы аргументов, т.е. если Base.specializations(@which f(...)) содержит специализации для рассматриваемого аргумента.

Break functions into multiple definitions

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

Вот пример "составной функции", которую на самом деле следует записать в виде нескольких определений:

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.

Write "type-stable" functions

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

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

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

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

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

Avoid changing the type of a variable

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

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

Separate kernel functions (aka, function barriers)

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

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, возникают при работе с данными неопределенного типа, например, с данными, загруженными из входного файла, которые могут содержать либо целые числа, либо числа с плавающей запятой, либо строки, либо что-то еще.

Types with values-as-parameters

Предположим, вы хотите создать 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; такой код будет очень медленным.

Теперь один очень хороший способ решения таких проблем — это использование function-barrier technique. Однако в некоторых случаях вы можете захотеть полностью устранить нестабильность типов. В таких случаях одним из подходов является передача размерности в качестве параметра, например через Val{T}() (см. "Value types"):

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

Юлия имеет специализированную версию 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)), либо уже указано в типовой области.

The dangers of abusing multiple dispatch (aka, more on types with values-as-parameters)

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

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", можно найти on the mailing list.

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

Access arrays in memory order, along columns

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

Pre-allocating outputs

Если ваша функция возвращает 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) для операций на месте с объединенными циклами и без временных массивов (см. dot syntax for vectorizing functions).

Use `MutableArithmetics` for more control over allocation for mutable arithmetic types

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

More dots: Fuse vectorized operations

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

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

Fewer dots: Unfuse certain intermediate broadcasts

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

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)

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

Consider using views for slices

В 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.

Copying data is not always bad

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

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

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

При условии, что достаточно памяти, стоимость копирования представления в массив компенсируется увеличением скорости от выполнения повторяющихся умножений матриц на непрерывном массиве.

Consider StaticArrays.jl for small fixed-size vector/matrix operations

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

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

Avoid string interpolation for I/O

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

println(file, "$a $b")

использовать:

println(file, a, " ", b)

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

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

против:

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

Optimize network I/O during parallel execution

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

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). 4d61726b646f776e2e436f64652822222c202266657463682229_40726566/4d61726b646f776e2e436f64652822222c2022776169742229_40726566 также выполняется последовательно, что приводит к общему ухудшению производительности.

Fix deprecation warnings

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

Tweaks

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

Избегайте ненужных массивов. Например, вместо 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).

Performance Annotations

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

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

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

Note

Хотя @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

Treat Subnormal Numbers as Zeros

Субнормальные числа, ранее называвшиеся denormal numbers, полезны во многих контекстах, но могут привести к снижению производительности на некотором оборудовании. Вызов 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). Самая важная характеристика 4d61726b646f776e2e436f64652822222c202240636f64655f7761726e747970652229_40726566 заключается в том, что неконкретные типы отображаются красным цветом; поскольку этот документ написан в 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 будут конкретными, и его производительность будет оптимальной. Однако есть обстоятельства, при которых такая эпhemeral типовая нестабильность может не иметь большого значения: например, если 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

Performance of captured variable

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

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.

Multithreading and linear algebra

Этот раздел относится к многопоточному коду 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. Однако это всего лишь общее правило, и лучший способ установить каждое количество потоков — это экспериментировать с вашим конкретным приложением.

Alternative linear algebra backends

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

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

Execution latency, package loading and package precompiling time

Reducing time to first plot etc.

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

Reducing package loading time

Снижение времени загрузки пакета обычно полезно. Общие рекомендации для разработчиков пакетов включают:

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

Инструмент @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.

Measure performance with @time and pay attention to memory allocation

Measure performance with `@time` and pay attention to memory allocation