Многомерные массивы

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

Массив — это коллекция объектов, хранящихся в многомерной сетке. В самом общем случае массив может содержать объекты типа Any. Для большинства вычислительных целей массивы должны содержать объекты более специфического типа, такие как Float64 или Int32.

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

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

Массивы

Основные функции

Создание и инициализация

Представлено множество функций для создания и инициализации массивов. В следующем списке таких функций вызовы с аргументом dims... могут принимать либо один кортеж размеров измерений, либо серию размеров измерений, переданных в качестве переменного числа аргументов. Большинство этих функций также принимают на вход первый аргумент T, который является типом элементов массива. Если тип T опущен, он будет по умолчанию Float64.

Синтаксис [A, B, C, ...] строит одномерный массив (вектор) из своих аргументов. Если все аргументы имеют общий тип повышения, они преобразуются к этому типу с использованием convert().

Конкатенация

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

Скалярные значения, переданные в эти функции, обрабатываются как массивы из 1 элемента.

Функции конкатенации используются так часто, что они имеют специальный синтаксис:

hvcat() выполняет конкатенацию по размерности 1 (с точками с запятой) и размерности 2 (с пробелами).

Инициализаторы массивов со строгим типом

Массив со строгим типом элементов может быть создан с помощью синтаксиса T[A, B, C, ...]. Это создаст одномерный массив с типом элементов T, инициализированный элементами A, B, C и т. д. Например, Any[x, y, z] создаёт гетерогенный массив, который может содержать любые значения.

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

julia> [[1 2] [3 4]]
1×4 Array{Int64,2}:
 1  2  3  4

julia> Int8[[1 2] [3 4]]
1×4 Array{Int8,2}:
 1  2  3  4

Понимания

Понимания обеспечивают общий и мощный способ создания массивов. Синтаксис понимания похож на обозначения создания множеств в математике:

A = [ F(x,y,...) for x=rx, y=ry, ... ]

Значение этого выражения заключается в том, что F(x,y,...) вычисляется с переменными x, y, и т. д., принимающими каждое значение из их заданного списка значений. Значения могут быть заданы как любой итерируемый объект, но обычно это будут диапазоны, такие как 1:n или 2:(n-1), или явные массивы значений, такие как [1.2, 3.4, 5.7]. Результатом является N-мерный плотный массив с измерениями, являющимися конкатенацией измерений диапазонов переменных rx, ry, и т. д., и каждое вычисление F(x,y,...) возвращает скаляр.

Следующий пример вычисляет взвешенное среднее значение текущего элемента и его левого и правого соседей вдоль одномерной сетки. :

julia> x = rand(8)
8-element Array{Float64,1}:
 0.843025
 0.869052
 0.365105
 0.699456
 0.977653
 0.994953
 0.41084
 0.809411

julia> [ 0.25*x[i-1] + 0.5*x[i] + 0.25*x[i+1] for i=2:length(x)-1 ]
6-element Array{Float64,1}:
 0.736559
 0.57468
 0.685417
 0.912429
 0.8446
 0.656511

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

Float32[ 0.25*x[i-1] + 0.5*x[i] + 0.25*x[i+1] for i=2:length(x)-1 ]

Генераторные выражения

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

julia> sum(1/n^2 for n=1:1000)
1.6439345666815615

При записи генераторного выражения с несколькими измерениями внутри списка аргументов необходимы круглые скобки для разделения генератора от последующих аргументов:

julia> map(tuple, 1/(i+j) for i=1:2, j=1:2, [1:4;])
ERROR: syntax: invalid iteration specification

Все выражения, разделённые запятыми, после for интерпретируются как диапазоны. Добавление круглых скобок позволяет нам добавить третий аргумент к map:

julia> map(tuple, (1/(i+j) for i=1:2, j=1:2), [1 3; 2 4])
2×2 Array{Tuple{Float64,Int64},2}:
 (0.5, 1)       (0.333333, 3)
 (0.333333, 2)  (0.25, 4)

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

julia> [(i,j) for i=1:3 for j=1:i]
6-element Array{Tuple{Int64,Int64},1}:
 (1, 1)
 (2, 1)
 (2, 2)
 (3, 1)
 (3, 2)
 (3, 3)

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

Сгенерированные значения можно отфильтровать с помощью ключевого слова if:

julia> [(i,j) for i=1:3 for j=1:i if i+j == 4]
2-element Array{Tuple{Int64,Int64},1}:
 (2, 2)
 (3, 1)

Индексирование

Общая синтаксис индексирования n-мерного массива A:

X = A[I_1, I_2, ..., I_n]

где каждый I_k может быть скалярным целым числом, массивом целых чисел или любым другим поддерживаемым индексом. Это включает Colon (:) для выбора всех индексов в рамках всего измерения, диапазоны вида a:c или a:b:c для выбора непрерывных или с шагом подсекций и массивы булевых значений для выбора элементов по своим true индексам.

Если все индексы являются скалярами, то результат X — это один элемент из массива A. В противном случае, X — это массив с тем же количеством измерений, что и сумма размерностей всех индексов.

Если все индексы являются векторами, например, то форма X будет (length(I_1), length(I_2), ..., length(I_n)), с местоположением (i_1, i_2, ..., i_n) из X содержащим значение A[I_1[i_1], I_2[i_2], ..., I_n[i_n]]. Если I_1 изменено на двумерную матрицу, то X становится n+1-мерным массивом формы (size(I_1, 1), size(I_1, 2), length(I_2), ..., length(I_n)). Матрица добавляет измерение. Местоположение (i_1, i_2, i_3, ..., i_{n+1}) содержит значение в A[I_1[i_1, i_2], I_2[i_3], ..., I_n[i_{n+1}]]. Все измерения, индексированные скалярами, отбрасываются. Например, результатом A[2, I, 3] является массив размером size(I). Его i-й элемент заполняется A[2, I[i], 3].

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

X = getindex(A, I_1, I_2, ..., I_n)

Пример:

julia> x = reshape(1:16, 4, 4)
4×4 Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> x[2:3, 2:end-1]
2×2 Array{Int64,2}:
 6  10
 7  11

julia> x[1, [2 3; 4 1]]
2×2 Array{Int64,2}:
  5  9
 13  1

Пустые диапазоны вида n:n-1 иногда используются для обозначения местоположения между индексами n-1 и n. Например, функция searchsorted() использует эту конвенцию для обозначения точки вставки значения, отсутствующего в отсортированном массиве:

julia> a = [1,2,5,6,7];

julia> searchsorted(a, 3)
3:2

Присваивание

Общая синтаксис присваивания значений в n-мерном массиве A:

A[I_1, I_2, ..., I_n] = X

где каждый I_k может быть скалярным целым числом, массивом целых чисел или любым другим поддерживаемым индексом. Это включает Colon (:) для выбора всех индексов в рамках всего измерения, диапазоны вида a:c или a:b:c для выбора непрерывных или с шагом подсекций и массивы булевых значений для выбора элементов по своим true индексам.

Если X — это массив, он должен иметь то же количество элементов, что и произведение длин индексов: prod(length(I_1), length(I_2), ..., length(I_n)). Значение в местоположении I_1[i_1], I_2[i_2], ..., I_n[i_n] из A перезаписывается значением X[i_1, i_2, ..., i_n]. Если X не является массивом, его значение записывается во все указанные местоположения A.

Точно так же, как и в индексировании, ключевое слово end может использоваться для представления последнего индекса каждого измерения в скобках индексирования, определяемого размером присваиваемого массива. Синтаксис индексированного присваивания без ключевого слова end эквивалентен вызову setindex!():

setindex!(A, X, I_1, I_2, ..., I_n)

Пример:

julia> x = collect(reshape(1:9, 3, 3))
3×3 Array{Int64,2}:
 1  4  7
 2  5  8
 3  6  9

julia> x[1:2, 2:3] = -1
-1

julia> x
3×3 Array{Int64,2}:
 1  -1  -1
 2  -1  -1
 3   6   9

Поддерживаемые типы индексов

В выражении A[I_1, I_2, ..., I_n], каждый I_k может быть скалярным индексом, массивом скалярных индексов или объектом, представляющим массив скалярных индексов и может быть преобразован в таковой функцией to_indices:

1. Скалярный индекс. По умолчанию это включает:

   • Целые числа, отличные от булевых

   • CartesianIndex{N}ы, которые ведут себя как N-кортеж целых чисел, охватывающих несколько измерений (см. ниже для более подробной информации)

2. Массив скалярных индексов. Это включает:

   • Векторы и многомерные массивы целых чисел

   • Пустые массивы, такие как [], которые не выбирают элементов

   • Rangeы вида a:c или a:b:c, которые выбирают непрерывные или с шагом подсекции от a до c (включительно)

   • Любой пользовательский массив скалярных индексов, являющийся подтипом AbstractArray

   • Массивы CartesianIndex{N} (см. ниже для более подробной информации)

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

   • Colon() (:), который представляет все индексы в рамках всего измерения или по всему массиву

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

Декартовы индексы

Специальный объект CartesianIndex{N} представляет скалярный индекс, который ведет себя как N-кортеж целых чисел, охватывающих несколько измерений. Например:

julia> A = reshape(1:32, 4, 4, 2);

julia> A[3, 2, 1]
7

julia> A[CartesianIndex(3, 2, 1)] == A[3, 2, 1] == 7
true

Взятый отдельно, это может показаться относительно тривиальным; CartesianIndex просто собирает несколько целых чисел в один объект, представляющий единственный многомерный индекс. Однако, в сочетании с другими формами индексирования и итераторами, которые возвращают CartesianIndexы, это может напрямую привести к очень элегантному и эффективному коду. См. Итерацию ниже, а для некоторых более продвинутых примеров см. эту запись блога о многомерных алгоритмах и итерации.

Также поддерживаются массивы CartesianIndex{N} . Они представляют собой набор скалярных индексов, каждый из которых охватывает N измерений, что позволяет использовать форму индексирования, иногда называемую точечным индексированием. Например, это позволяет получить доступ к диагональным элементам из первой "страницы" A из примера выше:

julia> page = A[:,:,1]
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> page[[CartesianIndex(1,1),
             CartesianIndex(2,2),
             CartesianIndex(3,3),
             CartesianIndex(4,4)]]
4-element Array{Int64,1}:
  1
  6
 11
 16

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

julia> A[CartesianIndex.(indices(A, 1), indices(A, 2)), 1]
4-element Array{Int64,1}:
  1
  6
 11
 16

julia> A[CartesianIndex.(indices(A, 1), indices(A, 2)), :]
4×2 Array{Int64,2}:
  1  17
  6  22
 11  27
 16  32

Логическое индексирование

Часто называемое логическим индексированием или индексированием с логической маской, индексирование булевым массивом выбирает элементы в тех индексах, где его значения равны true. Индексирование булевым вектором B по существу то же самое, что индексирование вектором целых чисел, который возвращается find(B). Аналогично, индексирование N-мерным булевым массивом по существу то же самое, что индексирование вектором CartesianIndex{N} с значениями true. Логический индекс должен быть вектором той же длины, что и измерение, в которое он индексируется, или он должен быть единственным предоставленным индексом и соответствовать размеру и размерности массива, в который он индексируется. Обычно более эффективно использовать булевы массивы в качестве индексов напрямую, вместо предварительного вызова find().

julia> x = reshape(1:16, 4, 4)
4×4 Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> x[[false, true, true, false], :]
2×4 Array{Int64,2}:
 2  6  10  14
 3  7  11  15

julia> mask = map(ispow2, x)
4×4 Array{Bool,2}:
  true  false  false  false
  true  false  false  false
 false  false  false  false
  true   true  false   true

julia> x[mask]
5-element Array{Int64,1}:
  1
  2
  4
  8
 16

Итерация

Рекомендуемые способы итерации по всему массиву

for a in A
    # Do something with the element a
end

for i in eachindex(A)
    # Do something with i and/or A[i]
end

Первая конструкция используется, когда вам нужно значение, но не индекс, каждого элемента. Во второй конструкции i будет Int, если A — это тип массива с быстрым линейным индексированием; в противном случае это будет CartesianIndex:

julia> A = rand(4,3);

julia> B = view(A, 1:3, 2:3);

julia> for i in eachindex(B)
           @show i
       end
i = CartesianIndex{2}((1, 1))
i = CartesianIndex{2}((2, 1))
i = CartesianIndex{2}((3, 1))
i = CartesianIndex{2}((1, 2))
i = CartesianIndex{2}((2, 2))
i = CartesianIndex{2}((3, 2))

В отличие от for i = 1:length(A), итерация с помощью eachindex обеспечивает эффективный способ итерации по любому типу массива.

Свойства массивов

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

Base.IndexStyle(::Type{<:MyArray}) = IndexLinear()

Это задаст использование целых чисел при итерации по MyArray. Если вы не укажете этот признак, будет использовано значение по умолчанию IndexCartesian().

Массивы и векторизованные операторы и функции

Для массивов поддерживаются следующие операторы:

1. Унарные арифметические – -, +

2. Бинарные арифметические – -, +, *, /, \, ^

3. Сравнения – ==, !=, ≈ (isapprox), ≉

Большинство бинарных арифметических операторов, перечисленных выше, также работают поэлементно, когда один из аргументов является скаляром: -, +, и * когда любой из аргументов является скаляром, а также / и \ когда знаменатель является скаляром. Например, [1, 2] + 3 == [4, 5] и [6, 4] / 2 == [3, 2].

Кроме того, для удобной векторизации математических и других операций, Julia предоставляет синтаксис с точкой f.(args...), например, sin.(x) или min.(x,y), для поэлементных операций над массивами или смешанными массивами и скалярами (операция распространения); эти операции обладают дополнительным преимуществом «слияния» в один цикл при комбинации с другими вызовами с точкой, например, sin.(cos.(x)).

Также каждый бинарный оператор поддерживает версию с точкой, которая может применяться к массивам (и комбинациям массивов и скаляров) в таких слитых операциях распространения, например, z .== sin.(x .* y).

Обратите внимание, что сравнения, такие как == , работают над целыми массивами, возвращая одно булево значение. Используйте операторы с точкой, такие как .== для поэлементных сравнений. (Для операций сравнения, таких как <, только поэлементная .< версия применима к массивам.)

Обратите также внимание на разницу между max.(a,b), который broadcast max() поэлементно по a и b, и maximum(a), который находит наибольшее значение в a. Такая же взаимосвязь существует для min.(a,b) и minimum(a).

Распространение

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

julia> a = rand(2,1); A = rand(2,3);

julia> repmat(a,1,3)+A
2×3 Array{Float64,2}:
 1.20813  1.82068  1.25387
 1.56851  1.86401  1.67846

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

julia> broadcast(+, a, A)
2×3 Array{Float64,2}:
 1.20813  1.82068  1.25387
 1.56851  1.86401  1.67846

julia> b = rand(1,2)
1×2 Array{Float64,2}:
 0.867535  0.00457906

julia> broadcast(+, a, b)
2×2 Array{Float64,2}:
 1.71056  0.847604
 1.73659  0.873631

Операторы с точкой, такие как .+ и .* эквивалентны вызовам broadcast (за исключением того, что они сливаются, как описано ниже). Также существует функция broadcast!() для указания явного назначения (к которому также можно получить доступ с помощью слияния присваиванием .= ), и функции broadcast_getindex() и broadcast_setindex!(), которые распространяют индексы перед индексированием. Кроме того, f.(args...) эквивалентно broadcast(f, args...), предоставляя удобный синтаксис для распространения любой функции (синтаксис с точкой). Вложенные «вызовы с точкой» f.(...) (включая вызовы .+ и т.д.) автоматически сливаются в один вызов broadcast.

Кроме того, broadcast() не ограничен массивами (см. документацию функции), он также обрабатывает кортежи и рассматривает любой аргумент, который не является массивом, кортежем или Ref (кроме Ptr) как «скаляр».

julia> convert.(Float32, [1, 2])
2-element Array{Float32,1}:
 1.0
 2.0

julia> ceil.((UInt8,), [1.2 3.4; 5.6 6.7])
2×2 Array{UInt8,2}:
 0x02  0x04
 0x06  0x07

julia> string.(1:3, ". ", ["First", "Second", "Third"])
3-element Array{String,1}:
 "1. First"
 "2. Second"
 "3. Third"

Реализация

Базовый тип массива в Julia – абстрактный тип AbstractArray{T,N}. Он параметризован количеством измерений N и типом элемента T. AbstractVector и AbstractMatrix – псевдонимы для случаев с 1 и 2 измерениями. Операции над объектами AbstractArray определены с помощью операторов и функций более высокого уровня, независимо от базовой реализации хранения. Эти операции в общем случае работают корректно в качестве обходного пути для любой конкретной реализации массива.

Тип AbstractArray включает в себя все объекты, напоминающие массив, и их реализации могут значительно отличаться от стандартных массивов. Например, элементы могут вычисляться по запросу, а не храниться. Однако любой конкретный тип AbstractArray{T,N} должен обычно реализовывать по крайней мере size(A) (возвращающий кортеж Int), getindex(A,i) и getindex(A,i1,...,iN); изменяемые массивы также должны реализовывать setindex!(). Рекомендуется, чтобы эти операции имели сложность практически постоянного времени, или, технически, сложность Õ(1), в противном случае некоторые функции массивов могут быть неожиданно медленными. Конкретные типы также обычно должны предоставлять метод similar(A,T=eltype(A),dims=size(A)), который используется для выделения аналогичного массива для copy() и других операций вне места. Независимо от того, как AbstractArray{T,N} представлен внутри, T – тип объекта, возвращаемого целочисленным индексированием (A[1, ..., 1], когда A не пусто), и N должна быть длиной кортежа, возвращенного size().

DenseArray – абстрактный подтип AbstractArray , предназначенный для включения всех массивов, расположенных в памяти с регулярными смещениями, которые, следовательно, могут быть переданы внешним функциям C и Fortran, ожидающим такое расположение памяти. Подтипы должны предоставлять метод stride(A,k), который возвращает «шаг» измерения k: увеличение индекса измерения k на 1 должно увеличить индекс i getindex(A,i) на stride(A,k). Если предоставляется метод преобразования указателя Base.unsafe_convert(Ptr{T}, A), расположение памяти должно соответствовать этим шагам аналогичным образом.

Тип Array – конкретный экземпляр DenseArray , где элементы хранятся в порядке по столбцам (см. дополнительные заметки в Рекомендациях по производительности). Vector и Matrix – псевдонимы для случаев с 1 и 2 измерениями. Специфические операции, такие как индексирование скаляром, присваивание и некоторые другие базовые операции, зависящие от хранения, – все, что нужно реализовать для Array, чтобы остальная часть библиотеки массивов могла быть реализована обобщенным способом.

SubArray – специализация AbstractArray , которая выполняет индексирование по ссылке, а не по копированию. SubArray создается с помощью функции view(), которая вызывается так же, как getindex() (с массивом и серией аргументов индекса). Результат view() выглядит так же, как результат getindex(), за исключением того, что данные остаются на месте. view() хранит входные векторы индексов в объекте SubArray , который позднее может быть использован для косвенного индексирования исходного массива. Разместив макрос @views перед выражением или блоком кода, любой срез array[...] в этом выражении будет преобразован для создания представления SubArray.

StridedVector и StridedMatrix – удобные псевдонимы, которые позволяют Julia вызывать более широкий спектр функций BLAS и LAPACK, передавая им либо Array , либо объекты SubArray , тем самым устраняя неэффективность из-за выделения памяти и копирования.

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

julia> a = rand(10,10)
10×10 Array{Float64,2}:
 0.561255   0.226678   0.203391  0.308912   …  0.750307  0.235023   0.217964
 0.718915   0.537192   0.556946  0.996234      0.666232  0.509423   0.660788
 0.493501   0.0565622  0.118392  0.493498      0.262048  0.940693   0.252965
 0.0470779  0.736979   0.264822  0.228787      0.161441  0.897023   0.567641
 0.343935   0.32327    0.795673  0.452242      0.468819  0.628507   0.511528
 0.935597   0.991511   0.571297  0.74485    …  0.84589   0.178834   0.284413
 0.160706   0.672252   0.133158  0.65554       0.371826  0.770628   0.0531208
 0.306617   0.836126   0.301198  0.0224702     0.39344   0.0370205  0.536062
 0.890947   0.168877   0.32002   0.486136      0.096078  0.172048   0.77672
 0.507762   0.573567   0.220124  0.165816      0.211049  0.433277   0.539476

julia> b = view(a, 2:2:8,2:2:4)
4×2 SubArray{Float64,2,Array{Float64,2},Tuple{StepRange{Int64,Int64},StepRange{Int64,Int64}},false}:
 0.537192  0.996234
 0.736979  0.228787
 0.991511  0.74485
 0.836126  0.0224702

julia> (q,r) = qr(b);

julia> q
4×2 Array{Float64,2}:
 -0.338809   0.78934
 -0.464815  -0.230274
 -0.625349   0.194538
 -0.527347  -0.534856

julia> r
2×2 Array{Float64,2}:
 -1.58553  -0.921517
  0.0       0.866567

Разреженные векторы и матрицы

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

Хранение разреженных матриц в формате Compressed Sparse Column (CSC)

В Julia разреженные матрицы хранятся в формате Compressed Sparse Column (CSC). Разреженные матрицы Julia имеют тип SparseMatrixCSC{Tv,Ti}, где Tv – тип хранимых значений, а Ti – целочисленный тип для хранения указателей столбцов и индексов строк. Внутреннее представление SparseMatrixCSC таково:

struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    m::Int                  # Number of rows
    n::Int                  # Number of columns
    colptr::Vector{Ti}      # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval::Vector{Ti}      # Row indices of stored values
    nzval::Vector{Tv}       # Stored values, typically nonzeros
end

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

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

Если у вас есть данные в формате CSC из другого приложения или библиотеки, и вы хотите импортировать их в Julia, убедитесь, что вы используете индексацию с началом с 1. Индексы строк в каждом столбце должны быть отсортированы. Если ваш SparseMatrixCSC объект содержит неотсортированные индексы строк, одним быстрым способом отсортировать их является двойной транспонирование.

В некоторых приложениях удобно хранить явные нулевые значения в SparseMatrixCSC. Эти значения приняты функциями в Base (но нет гарантии, что они будут сохранены в изменяющих операциях). Такие явно хранящиеся нули рассматриваются как структурные ненулевые элементы многими процедурами. Функция nnz() возвращает количество элементов, явно хранящихся в структуре разреженных данных, включая структурные ненулевые элементы. Для подсчёта точного количества числовых ненулевых элементов используйте countnz(), которая проверяет каждый хранящийся элемент разреженной матрицы. dropzeros() и встроенная dropzeros!() могут использоваться для удаления хранящихся нулей из разреженной матрицы.

julia> A = sparse([1, 2, 3], [1, 2, 3], [0, 2, 0])
3×3 SparseMatrixCSC{Int64,Int64} with 3 stored entries:
  [1, 1]  =  0
  [2, 2]  =  2
  [3, 3]  =  0

julia> dropzeros(A)
3×3 SparseMatrixCSC{Int64,Int64} with 1 stored entry:
  [2, 2]  =  2

Хранение разреженных векторов

Разреженные векторы хранятся в аналоге сжатого разреженного столбца для разреженных матриц. В Julia разреженные векторы имеют тип SparseVector{Tv,Ti}, где Tv — тип хранимых значений, а Ti — целочисленный тип для индексов. Внутреннее представление таково:

struct SparseVector{Tv,Ti<:Integer} <: AbstractSparseVector{Tv,Ti}
    n::Int              # Length of the sparse vector
    nzind::Vector{Ti}   # Indices of stored values
    nzval::Vector{Tv}   # Stored values, typically nonzeros
end

Как и для SparseMatrixCSC, тип SparseVector также может содержать явно хранящиеся нули. (См. Хранение разреженных матриц.)

Конструкторы разреженных векторов и матриц

Самый простой способ создания разреженных массивов — использование функций, эквивалентных функциям zeros() и eye(), которые Julia предоставляет для работы с плотно заполненными массивами. Для создания разреженных массивов вместо этого можно использовать те же имена с префиксом sp:

julia> spzeros(3)
3-element SparseVector{Float64,Int64} with 0 stored entries

julia> speye(3,5)
3×5 SparseMatrixCSC{Float64,Int64} with 3 stored entries:
  [1, 1]  =  1.0
  [2, 2]  =  1.0
  [3, 3]  =  1.0

Функция sparse() часто является удобным способом построения разреженных массивов. Например, для построения разреженной матрицы можно ввести вектор I индексов строк, вектор J индексов столбцов и вектор V хранимых значений (это также известно как формат COO (координат)). sparse(I,J,V) затем строит разреженную матрицу таким образом, что S[I[k], J[k]] = V[k]. Эквивалентный конструктор разреженного вектора — sparsevec, который принимает вектор индексов (строк) I и вектор V с хранимыми значениями и строит разреженный вектор R такой, что R[I[k]] = V[k].

julia> I = [1, 4, 3, 5]; J = [4, 7, 18, 9]; V = [1, 2, -5, 3];

julia> S = sparse(I,J,V)
5×18 SparseMatrixCSC{Int64,Int64} with 4 stored entries:
  [1 ,  4]  =  1
  [4 ,  7]  =  2
  [5 ,  9]  =  3
  [3 , 18]  =  -5

julia> R = sparsevec(I,V)
5-element SparseVector{Int64,Int64} with 4 stored entries:
  [1]  =  1
  [3]  =  -5
  [4]  =  2
  [5]  =  3

Обратной операцией к функциям sparse() и sparsevec является findnz(), которая извлекает входные данные, используемые для создания разреженного массива. Также есть функция findn, которая возвращает только векторы индексов.

julia> findnz(S)
([1, 4, 5, 3], [4, 7, 9, 18], [1, 2, 3, -5])

julia> findn(S)
([1, 4, 5, 3], [4, 7, 9, 18])

julia> findnz(R)
([1, 3, 4, 5], [1, -5, 2, 3])

julia> findn(R)
4-element Array{Int64,1}:
 1
 3
 4
 5

Другой способ создать разреженный массив — преобразовать плотно заполненный массив в разреженный массив с помощью функции sparse():

julia> sparse(eye(5))
5×5 SparseMatrixCSC{Float64,Int64} with 5 stored entries:
  [1, 1]  =  1.0
  [2, 2]  =  1.0
  [3, 3]  =  1.0
  [4, 4]  =  1.0
  [5, 5]  =  1.0

julia> sparse([1.0, 0.0, 1.0])
3-element SparseVector{Float64,Int64} with 2 stored entries:
  [1]  =  1.0
  [3]  =  1.0

Вы можете перейти в обратном направлении, используя конструктор Array. Функция issparse() может использоваться для проверки, является ли матрица разреженной.

julia> issparse(speye(5))
true

Операции с разреженными матрицами

Арифметические операции с разреженными матрицами также работают так же, как и с плотно заполненными матрицами. Индексирование, присваивание и конкатенация разреженных матриц работают так же, как и с плотно заполненными матрицами. Операции индексирования, особенно присваивание, являются дорогостоящими, когда выполняются по одному элементу за раз. Во многих случаях может быть лучше преобразовать разреженную матрицу в формат (I,J,V) с помощью findnz(), обработать значения или структуру в плотно заполненных векторах (I,J,V), а затем восстановить разреженную матрицу.

Соответствие плотно заполненных и разреженных методов

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

Подробности можно найти в разделе «Разреженные векторы и матрицы» справочника стандартной библиотеки.