Julia: эффективный и простой в использовании язык программирования для числовых вычислений/оптимизации.

Сегодня я отдыхаю от своего напряженного графика и кратко расскажу о Джулии, очень простом в использовании языке программирования. Лично для тех, кто занимается оптимизацией/числовыми расчетами, это лучший язык для использования на данный момент. Что касается меня, я не использовал Matlab с тех пор, как использовал Julia в течение 3 лет, и я думаю, что Python в основном бесполезен.

The Julia Language

1. Введение в числовые вычисления в Джулии

Во-первых, основанный на LLVM JIT-компилятор Джулии обеспечивает менее сложную среду записи (по сравнению с C) и по-прежнему обеспечивает производительность, аналогичную C по времени числовых вычислений.Julia Micro-Benchmarks). Протестированные операции включают $\pi$ Суммирование рядов, рекурсивное суммирование/быстрая сортировка Фибоначчи, вывод в файл, умножение матриц/собственные значения и т. д.

Пример 1: время достижения для простого случайного блуждания

Функция walk(L) определяет x=2 является отправной точкой, а правая граница x=L , левая граница x=1 Простое случайное блуждание . В этом примере мы в основном рассматриваем время достижения левой границы.Функция times(L,N) возвращает время, когда walk(L) впервые достигает левой границы.Повторение(L,N,M) функция позволяет нам повторять прогулку несколько раз (т. е. симуляция Монте-Карло) и возвращать среднее время попадания. В частности, всего имеется M симуляций, и каждая симуляция повторяется N раз и усредняется.

function walk(L=50)
    t = 2
    x = 2
    while x != 1
        if rand() < 0.5
            x += 1
        else
            x -= 1
        end 
        if x > L
            x = L
        end
        t += 1
    end
    return t
end

function run(L, N)
    times = Int64[]
    for i in 1:N
        push!(times, walk(L))
    end
    return times
end

function repeat(L, N, M)
    [mean(run(L, N)) for i in 1:M]
end

Среднее время попадания, возвращенное нашей симуляцией Монте-Карло, составило 99,90814250000001, и была выведена соответствующая гистограмма pmf.

data = repeat(50, 1000, 10000);
using Plots
mean(data)
histogram(data)

Тогда мы можем видеть, что pmf больше похож на распределение Гаусса, а среднее сосредоточено во. Конечно, это очевидно, потому что мы можем вычислить это вручную для проверки. определение $\mu_i=\mathbb{E}[\min\{ n\geq 0| X_n = 1 \}|X_0 = i ]$ , так легко пройти $\mu_1=2,\mu_2=1+\frac{1}{2}\mu_1+\frac{1}{2}\mu_3,\mu_3=1+\frac{1}{2}\mu_2+\frac{1}{2}\mu_4,\ldots,\mu_{L-1}=1+\frac{1}{2}\mu_{L-2}+\frac{1}{2}\mu_L,\mu_L=1+\frac{1}{2}\mu_{L-1}+\frac{1}{2}\mu_L$ получить $\mu_2=2L.$ Поскольку распределение выглядит как гауссово, это из-за центральной предельной теоремы.

Поскольку это случайное блуждание также можно рассматривать как цепь Маркова, мы также можем рассчитать его стационарное распределение через произведение матрицы перехода состояния, Ну, это в основном для того, чтобы показать больше команд Джулии и сравнить понимание массива и (разреженную) скорость умножение матриц. Мы используем пакет BenchmarkTools для сравнения скорости этих трех алгоритмов. позволять L=500 и повторить 1000 раз (tmax=1000).

# Array comprehensions
function evolve(p)
    L = length(p)
    
    p_new = vec(zeros(L,1))
    p_new[1] = p[1] + p[2]/2
    p_new[2] = p[3] / 2
    
    for i in 3:L-1
        p_new[i] = (p[i-1] + p[i+1]) / 2
    end
    
    p_new[L] = (p[L-1] + p[L]) / 2
    
    return p_new
end

@benchmark begin
    p = [i == 2 ? 1//1 : 0//1 for i in 1:L];
    P = zeros(t_max,L)
    P[1,:] = evolve(p)
    for i = 2:t_max
        P[i,:] = evolve(P[i-1,:])
    end
end

# Dense matrix
@benchmark begin
    p = [i == 2 ? 1//1 : 0//1 for i in 1:L];
    P = zeros(t_max,L)
    Trans = zeros(L,L)
    Trans[1,1] = 1
    Trans[1,2] = 0.5
    Trans[2,3] = 0.5
    for i = 3:L-1
        Trans[i,i-1] = 0.5
        Trans[i,i+1] = 0.5
    end
    Trans[L,L-1] = 0.5
    Trans[L,L] = 0.5
    P[1,:] = Trans * vec(p)
    for i = 2:t_max
        P[i,:] = Trans * vec(P[i-1,:])
    end
end

# Sparse matrix
@benchmark begin
    t_max = 1000
    p = [i == 2 ? 1//1 : 0//1 for i in 1:L];
    P = zeros(t_max,L)
    dl = 0.5 * vec(ones(L-1,1))
    dl[1] = 0
    d = vec(zeros(L,1))
    d[1] = 1
    d[L] = 0.5
    du = 0.5 * vec(ones(L-1,1))
    Trans = Tridiagonal(dl,d,du)
    P[1,:] = Trans * vec(p)
    for i = 2:t_max
        P[i,:] = Trans * vec(P[i-1,:])
    end
end

Обратите внимание, что в разреженной матрице мы используем трехдиагональную матрицу (Tridiagonal) структуру данных, и тест-инструменты сделали сотни выборок, а затем вывели различную статистику времени работы блока кода. Затем мы знаем, что разреженные матрицы быстрее, чем понимание массива быстрее, чем наивные плотные матрицы. В численных вычислениях всегда важно использовать наиболее эффективные структуры данных и кодовые команды. Что касается того, если вы хотите проверить (проверить работоспособность), дают ли эти три блока кода согласованные результаты вывода, вы можете использовать команду @test_ приблизительно_eq, которая не будет повторяться здесь.

Пример 2: Разреженная матрица, пользовательская структура данных, распределенная структура данных (предварительно)

В приведенном выше примере мы видели применение специального класса разреженных матриц — типа трехдиагональной матрицы. Здесь мы указываем, что более общий тип разреженной матрицы называется sparse() в Джулии. Простой пример выглядит следующим образом:

G = sparse(vec([1;2;3;1;1;2;3]),vec([1;2;3;2;3;1;1]),vec([1.;3;5;2;2;2;2]))

Мы видим, что этот тип матрицы будет хранить только ненулевые элементы и соответствующие им значения, что очень похоже на Matlab. Затем мы пытаемся настроить новый тип матрицы (структуру данных), мы называем их матрицей SymArrow, своего рода симметричной матрицей, за исключением диагонали и первой строки/столбца, другие элементы являются разреженной матрицей 0.

type SymArrowFloat
    Diag::Array{Float64,1}
    Row::Array{Float64,1}
end

import Base: full 
Base.full(A::SymArrowFloat) = full(sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ))

import Base: promote_rule, convert, show
Base.show(io::IO, A::SymArrowFloat)=print(io,sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ) )

import Base.+
+(A::SymArrowFloat,B::SymArrowFloat) = SymArrowFloat(A.Diag+B.Diag, A.Row+B.Row)

import Base.*
*(A::SymArrowFloat,v::Array{Float64,1}) = [dot(vec([A.Diag[1];A.Row]),vec(v));[v[1]*A.Row[i-1]+v[i]*A.Diag[i] for i=2:size(vec(v),1)]]

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

# Function full works
G = SymArrowFloat([1.;3;5],[2;2.])
full(G)
# Function show works
G
# + works
G+G
# * works
G*rand(3)

Результаты работы здесь не выкладываются, можете сами попробовать и все работает. Здесь мы отмечаем, что если вы хотите указать рациональные числа, комплексные числа могут быть определены отдельно (на самом деле, это в основном просто для отображения числового поля Джулии типа 233).

type SymArrow{T}
    Diag::Array{T,1}
    Row::Array{T,1} 
end

import Base: full 
Base.full{T}(A::SymArrow{T}) = full(sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ))

import Base.+
+{T}(A::SymArrow{T},B::SymArrow{T}) = SymArrow{T}(A.Diag+B.Diag, A.Row+B.Row)

import Base.*
*{T}(A::SymArrow{T},v::Array{T,1}) = [dot(vec([A.Diag[1];A.Row]),vec(v));[v[1]*A.Row[i-1]+v[i]*A.Diag[i] for i=2:size(vec(v),1)]]

import Base: promote_rule, convert, show
Base.show{T}(io::IO, A::SymArrow{T})=print(io,sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ) )

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

G = SymArrow{Rational}([1.;3;5],[2;2.])

G = SymArrow{Complex}([1.+0im;3;5],[2;2.-0im])

G = SymArrow{BigFloat}([1.;3;5],[2;2.])
G*[BigFloat(1);2;3]

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

type SymArrow2{T} <: AbstractMatrix{T}
    Diag::Array{T,1}
    Row::Array{T,1}  
end

import Base: size, getindex

size{T}(A::SymArrow2{T}) = (length(A.Diag), length(A.Diag))

function getindex{T}(A::SymArrow2{T}, i, j)  
    if i == j
        return A.Diag[i]
    elseif i == 1
        return A.Row[j-1]
    elseif j == 1
        return A.Row[i-1]
    else
        return zero(T)  # otherwise return zero of type T
    end
end

В конце этого примера давайте рассмотрим распределенные структуры данных, поддерживаемые Julia (здесь кратко представлены только массивы), и некоторые связанные с ними распределенные операции. Что же касается реального использования Julia для параллельных вычислений, то в этой статье мы его обсуждать не будем, возможно, в будущем он станет более популярным. Здесь мы рассматриваем только пакет DistributedArrays.jl. Примером является распределенная кумулятивная сумма по массиву. Мы используем 4 ядра для распределенного выполнения операций.

using DistributedArrays
nprocs()  ==  1  &&  addprocs(4)
xd = @DArray [i for i = 0:39]
[DistributedArrays.chunk(xd, i) for i = 1:4]

Мы видим, что массив автоматически разбивается на 4 части и распределяется по соответствующим областям. Чтобы выполнять параллельные операции, нам нужно определить функцию prefix(), в основном часть map в структуре map-reduce.

@everywhere  function prefix!(op, v, v0 = 0)
    v[1] += v0
    for i = 2:length(v)
        v[i] = op(v[i-1], v[i])
    end
    return v
end

@everywhere function prefixlag!(op, v, v0 = 0)
    vi = v[1]
    v[1] = v0
    for i = 2:length(v)
        tmp  = op(v[i - 1], vi)
        vi   = v[i]
        v[i] = tmp
    end
    return v
end

Затем мы можем получить cumsum распределенно.

yd1 = map_localparts(t -> [sum(t)], xd)
yd2 = prefixlag!(+, Array(yd1))
yd3 = map_localparts((t,s) -> prefix!(+, t, s[1]), xd, distribute(yd2))

Пример 3: Разложение по единственному значению (SVD)

Этот пример является исходным кодом в моем ответе.Цинь Ханьчжан: Где отражается уменьшение размерности СВД?

using Images, Colors, ImageCore, ImageView
img = load("Albert Einstein-rare-pics56.jpg")
img_RGB = float(Array(channelview(img)))
R_pixel = img_RGB[1,:,:];
G_pixel = img_RGB[2,:,:];
B_pixel = img_RGB[3,:,:];
U1,S1,V1 = svd(R_pixel);
U2,S2,V2 = svd(G_pixel);
U3,S3,V3 = svd(B_pixel);
U1[:,1:100]*diagm(S1[1:100])*V1[1:100,:]

# Full SVD: perfect reconstruction
Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))
# 25% SVD reconstruction: near perfect
S1[101:400] .= 0
S2[101:400] .= 0
S3[101:400] .= 0
Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))

# 2.5% SVD reconstruction: obviously blurred
S1[11:400] .= 0
S2[11:400] .= 0
S3[11:400] .= 0
Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))

# 1% SVD reconstruction: heavily blurred
S1[5:400] .= 0
S2[5:400] .= 0
S3[5:400] .= 0
Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))

# only 2 nonzero entries: stripes
S1[3:400] .= 0
S2[3:400] .= 0
S3[3:400] .= 0
Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))

2. Введение в JuMP: линейное программирование, конусное программирование, смешанное целочисленное программирование.

Помимо численных расчетов, язык Julia очень удобен для оптимизаторов, аналог ЯЛМИП в Matlab (yalmip.github.io/), предоставляет входной интерфейс модели оптимизации высокого уровня и вызывает различные решатели для решения модели оптимизации. В Юлии это в основном достигается за счет пакета JuMP.

Список поддерживаемых в настоящее время решателей: (см.Optimization packages for the Julia language)

С точки зрения скорости моделирования он строго лучше, чем другие интерфейсы AML (язык алгебраического моделирования) на рынке, и не уступает коммерческому AML!

Источник: Даннинг, Иэн, Джоуи Хьючетт и Майлз Любин, «JuMP: язык моделирования для математической оптимизации», SIAM Review 59.2 (2017): 295–320.

Пример 1: Машина опорных векторов (SVM) [линейное программирование]

В этом примере я даю реализацию SVM в Julia с использованием линейного программирования, линейного ядра и ядра RBF с двумя функциями ядра. В частности, здесь рассматривается двойное представление SVM с мягкой маржой, а соответствующая стоимость резервной переменной равна

using PyPlot, JuMP

function linear_K(X)
    K = zeros(size(X,1),size(X,1))
    for i = 1:size(X,1)
        for j = 1:size(X,1)
            K[i,j] = dot(vec(X[i,:]), vec(X[j,:]))
        end
    end
    return K
end

function RBF_K(X,γ)
    K = zeros(size(X,1),size(X,1))
    for i = 1:size(X,1)
        for j = 1:size(X,1)
            K[i,j] = exp( -γ*sum(( X[i,:]-X[j,:] ).^2) )
        end
    end
    return K
end

function SVM_K(X,Y,C,K_type,γ=1)
    n = length(Y)
    M_dual = Model(solver = GurobiSolver(OutputFlag=0))
    @variable(M_dual, 0<=α[1:n]<=C)
    @constraint(M_dual, sum(Y.*α)==0)
    if K_type == "linear"
        @objective(M_dual, Max, -0.5*sum(α'*linear_K(X)*α) + sum(α))
    elseif K_type == "RBF"
        @objective(M_dual, Max, -0.5*sum(α'*RBF_K(X,γ)*α) + sum(α))
    else
        error("Wrong kernel type!")
    end
    solve(M_dual)
    # return optimal α and margin
    return getvalue(α), 1/(getobjectivevalue(M_dual)-sum(getvalue(α)))*(-2)
end

Обратите внимание, что в модели JuMP здесь я вызываю Gurobi в качестве решателя.

Пример 2: Расчет стабильного числа графа [полуположительно определенное программирование]

Для неориентированного графа G = (V, E) стабильным множеством называется подмножество узлов V, где никакие два узла не соединены. Стабильное число относится к максимальному количеству узлов во всех стабильных множествах на G. Эту проблему можно свести к полуопределенной задаче программирования:

$\begin{align} \alpha(G) = & \max_X~ \text{tr}(ee^TX) \\ &\text{s.t. } X_{ij} = 0,~ \text{if } (i,j)\in E,\\ & \text{tr}(X) = 1,\\ & X \succeq 0. \end{align}$

Код Джулии выглядит следующим образом. Суть в том, чтобы использовать макрос @SDconstraint, и обратите внимание, что здесь я вызываю Mosek в качестве решателя. Вот несколько конкретных примеров:Введение в копозитивное программирование

using JuMP, Mosek
function StableSDP(A)
    n = size(A,1)
    StableDD = Model(solver = MosekSolver(MSK_IPAR_LOG = 0))
    @variable(StableDD, X[1:n,1:n])
    @SDconstraint(StableDD, X>=0)
    @constraint(StableDD, X.>=0)
    @constraint(StableDD, sum(sum(A.*X)) == 0 )
    @constraint(StableDD, sum(X[i,i] for i=1:n) == 1)
    @objective(StableDD, Max, sum(sum(X)))
    solve(StableDD)
    return getobjectivevalue(StableDD)
end

Пример 3: Gurobi Advanced [Смешанное целочисленное программирование]

Здесь мы отмечаем, что, используя интерфейс JuMP, мы по-прежнему имеем возможность относительно свободно взаимодействовать с Gurobi. Например, мы можем добавить ленивое ограничение (то есть указать только часть ограничения в начале, а остальные ограничения добавлять постепенно по мере решения) или пользовательское ограничение, которое можно перемежать в процессе ветвления и привязки. И Gurobi также поддерживает теплый старт, то есть пользователь может указать начальное значение (например, полученное с помощью некоторой эвристики, ускоряющей ветвление и привязку); подробные руководства по этому поводу можно увидетьSolver Callbacks.

Точно так же Gurobi также можно использовать для генерации столбцов (аналогично ленивому ограничению, за исключением того, что в начале добавляются только некоторые переменные, а другие переменные добавляются постепенно по мере выполнения решения). К сожалению, Gurobi еще не поддерживает ветвь и цену (то есть генерацию столбцов в целочисленном программировании), поэтому в настоящее время CG можно использовать только в задачах непрерывной оптимизации. В этом отношении,Chiwei Yan- JuliaTutorialУчебное пособие по проблеме раскроя заготовки является отличным учебным материалом.

3. Интерактивное программирование

Отметим, что, подобно IPython, пакет IJulia в Julia позволяет выполнять всю компиляцию Julia в блокноте Jupyter. Лично мне нравится программировать Джулию в этой среде (см.:JuliaLang/IJulia.jl). Конечно, я знаю, что многим нравится Juno's:Juno, the Interactive Development EnvironmentМожет быть, этот больше похож на кодового фермера, ххх

Далее у нас может быть множество интерактивных операций, которые в основном реализованы в Джулии через пакет Interact. Например, мы можем настроить ползунки, кнопки и т. д. для взаимодействия с параметрической кривой. Посмотрите видео ниже для примера илиJuliaGizmos/Interact.jl.

Interact in Julia

Что касается всех видов графиков, я предпочитаю использовать пакет Plots. Начало работы можно найти по адресу:Plots - powerful convenience for visualization in Julia.

В-четвертых, напишите в конце

Естественно, в этой статье дано лишь несколько примеров и базовое введение в язык программирования Julia. Если вы просто хотите, чтобы удобный язык вызывал решатель, или в основном выполнял численные вычисления, или более продвинутый эксперт по алгоритмам оптимизации, больше практикуйтесь и учитесь, пока «пачкаете руки». Я думаю, что это всегда наиболее эффективно.

Во-первых, это большое количество учебных ресурсов, представленных на официальном сайте Юлии: включая видео, конкретные примеры и различные учебные пособия в формате pdf.

Learning Julia

Личная домашняя страница курса численных вычислений профессора Эдельмана, одного из соучредителей Julia, также является очень хорошей ссылкой: (в основном ссылка на Github)

Modern Numerical Computing with Julia

Наконец, немного личной надежды. Студенты, знакомые с Джулией, также могут публиковать соответствующие руководства на Zhihu, такие как некоторые более продвинутые приложения (надежная оптимизация, глубокое обучение, символьные вычисления, крупномасштабные распределенные вычисления и т. д.), мою личную колонку и«Оперативное планирование ИЛИ стратегия» Операционные исследования в эпоху искусственного интеллекта больших данныхКолонки также могут способствовать продвижению публикации. Я надеюсь, что Джулия и сообщество операционных исследований/оптимизации/ИЛИ станут сильнее~

Связанные ответы Quora:

Какие болевые точки Джулия решает для C++/Python/Matlab? Что вы думаете о новом языке Julia? Каковы его перспективы?