Spec-Zone.ru › Julia 0.5

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

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

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

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

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

Массивы

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

Функция Описание
eltype(A) тип элементов, содержащихся в A
length(A) количество элементов в A
ndims(A) количество измерений A
size(A) кортеж, содержащий размеры A
size(A,n) размер A вдоль определенного измерения
indices(A) кортеж, содержащий допустимые индексы A
indices(A,n) диапазон, выражающий допустимые индексы вдоль измерения n
eachindex(A) эффективный итератор для посещения каждой позиции в A
stride(A,k) шаг (линейное расстояние между смежными элементами) вдоль измерения k
strides(A) кортеж шагов в каждом измерении

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

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

Функция Описание
Array{type}(dims...) неинициализированный плотный массив
zeros(type, dims...) массив всех нулей заданного типа, по умолчанию Float64 если type не указан
zeros(A) массив всех нулей того же типа элементов и формы, что и A
ones(type, dims...) массив всех единиц заданного типа, по умолчанию Float64 если type не указан
ones(A) массив всех единиц того же типа элементов и формы, что и A
trues(dims...) массив Bool со всеми значениями true
trues(A) массив Bool со всеми значениями true и формой A
falses(dims...) массив Bool со всеми значениями false
falses(A) массив Bool со всеми значениями false и формой A
reshape(A, dims...) массив с теми же данными, что и данный массив, но с другими размерами.
copy(A) копия A
deepcopy(A) копия A, рекурсивно копирующая её элементы
similar(A, element_type, dims...) неинициализированный массив того же типа, что и заданный массив (плотный, разреженный и т.д.), но со специфицированным типом элементов и размерами. Второй и третий аргументы оба необязательны, по умолчанию принимая тип элемента и размеры A в случае опускания.
reinterpret(type, A) массив с теми же двоичными данными, что и заданный массив, но с заданным типом элемента
rand(dims) Array из Float64 с случайными, независимыми и одинаково распределёнными значениями в полуоткрытом интервале \([0, 1)\)
randn(dims) Array из Float64 со случайными, независимыми и стандартно нормально распределёнными случайными значениями
eye(n) матрица n на n единичная матрица
eye(m, n) матрица m на n единичная матрица
linspace(start, stop, n) диапазон из n линейно размещённых элементов от start до stop
fill!(A, x) заполнить массив A значением x
fill(x, dims) создать массив, заполненный значением x
[1] iid, независимые и одинаково распределённые.

Синтаксис [A, B, C, ...] создаёт одномерный массив (вектор) из своих аргументов.

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

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

Функция Описание
cat(k, A...) склеивание входных n-мерных массивов вдоль размерности k
vcat(A...) краткая запись для cat(1, A...)
hcat(A...) краткая запись для cat(2, A...)

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

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

Выражение Вызовы
[A; B; C; ...] vcat()
[A B C ...] hcat()
[A B; C D; ...] hvcat()

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 может быть:

  1. Целое скалярное значение
  2. Диапазон в форме a:b, или a:b:c
  3. Диапазон или Colon() для выбора целых размерностей
  4. Произвольный целочисленный массив, включая пустой массив []
  5. Булевый массив для выбора вектора элементов в его 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].

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

Кроме того, отдельные элементы многомерного массива можно индексировать как x = A[I], где I является CartesianIndex. Это эффективно ведёт себя как n-кортеж целых чисел, охватывающий несколько размерностей A. См. Итерацию ниже.

В качестве особой части этого синтаксиса ключевое слово 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[map(ispow2, x)]
5-element Array{Int64,1}:
  1
  2
  4
  8
 16

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 может быть:

  1. Целое скалярное значение
  2. Диапазон в форме a:b, или a:b:c
  3. Диапазон или Colon() для выбора целых размерностей
  4. Произвольный целочисленный массив, включая пустой массив []
  5. Булевый массив для выбора элементов в его 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.

Булевый массив, используемый в качестве индекса, ведёт себя как в getindex(), как будто он сначала преобразуется с помощью find().

Синтаксис присваивания индексов эквивалентен вызову 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

Итерация

Рекомендованные способы перебора всего массива:

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:

A = rand(4,3)
B = view(A, 1:3, 2:3)
julia> for i in eachindex(B)
           @show i
       end
       i = Base.IteratorsMD.CartesianIndex_2(1,1)
       i = Base.IteratorsMD.CartesianIndex_2(2,1)
       i = Base.IteratorsMD.CartesianIndex_2(3,1)
       i = Base.IteratorsMD.CartesianIndex_2(1,2)
       i = Base.IteratorsMD.CartesianIndex_2(2,2)
       i = Base.IteratorsMD.CartesianIndex_2(3,2)

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

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

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

Base.linearindexing{T<:MyArray}(::Type{T}) = LinearFast()

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

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

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

  1. Унарная арифметика — -, +, !
  2. Бинарная арифметика — +, -, *, .*, /, ./, \, .\, ^, .^, div, mod
  3. Сравнение — .==, .!=, .<, .<=, .>, .>=
  4. Унарная булева или побитовая — ~
  5. Бинарная булева или побитовая — &, |, $

Некоторые операторы без точек всё равно работают поэлементно, когда один аргумент является скаляром. Эти операторы — *, +, -, и побитовые операторы. Операторы / и \ работают поэлементно, когда знаменатель является скаляром.

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

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

abs abs2 angle cbrt
airy airyai airyaiprime airybi airybiprime airyprime
acos acosh asin asinh atan atan2 atanh
acsc acsch asec asech acot acoth
cos  cospi cosh  sin  sinpi sinh  tan  tanh  sinc  cosc
csc  csch  sec  sech  cot  coth
acosd asind atand asecd acscd acotd
cosd  sind  tand  secd  cscd  cotd
besselh besseli besselj besselj0 besselj1 besselk bessely bessely0 bessely1
exp  erf  erfc  erfinv erfcinv exp2  expm1
beta dawson digamma erfcx erfi
exponent eta zeta gamma
hankelh1 hankelh2
 ceil  floor  round  trunc
isfinite isinf isnan
lbeta lfact lgamma
log log10 log1p log2
copysign max min significand
sqrt hypot

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

Julia предоставляет макросы @vectorize_1arg() и @vectorize_2arg() для автоматической векторизации любой функции от одного или двух аргументов соответственно. Каждый из них принимает два аргумента: Type аргумента (обычно выбирается максимально общим) и имя векторизуемой функции. Вот простой пример:

julia> square(x) = x^2
square (generic function with 1 method)

julia> @vectorize_1arg Number square
square (generic function with 2 methods)

julia> methods(square)
# 2 methods for generic function "square":
square{T<:Number}(x::AbstractArray{T,N<:Any}) at operators.jl:555
square(x) at none:1

julia> square([1 2 4; 5 6 7])
2×3 Array{Int64,2}:
  1   4  16
 25  36  49

Векторизация

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

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_getindex() и broadcast_setindex!(), которые векторизуют индексы перед индексированием. Кроме того, f.(args...) эквивалентно broadcast(f, args...), обеспечивая удобный синтаксис для векторизации любой функции (Точечный синтаксис для векторизации функций).

Реализация

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

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

Разряженные матрицы

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

Хранение в формате Compressed Sparse Column (CSC)

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

type 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 values of nonzeros
    nzval::Vector{Tv}       # Nonzero values
end

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

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

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

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

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

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

julia> spzeros(3,5)
3×5 sparse matrix with 0 Float64 nonzero entries

julia> speye(3,5)
3×5 sparse matrix with 3 Float64 nonzero entries:
        [1, 1]  =  1.0
        [2, 2]  =  1.0
        [3, 3]  =  1.0

Функция sparse() часто является удобным способом построения разреженных матриц. Она принимает на вход вектор I индексов строк, вектор J индексов столбцов и вектор V ненулевых значений. sparse(I,J,V) создаёт разреженную матрицу, такую что S[I[k], J[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 sparse matrix with 4 Int64 nonzero entries:
        [1 ,  4]  =  1
        [4 ,  7]  =  2
        [5 ,  9]  =  3
        [3 , 18]  =  -5

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

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

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

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

julia> sparse(eye(5))
5×5 sparse matrix with 5 Float64 nonzero entries:
        [1, 1]  =  1.0
        [2, 2]  =  1.0
        [3, 3]  =  1.0
        [4, 4]  =  1.0
        [5, 5]  =  1.0

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

julia> issparse(speye(5))
true

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

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

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

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

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

Разреженная Плотная Описание
spzeros(m,n) zeros(m,n) Создаёт m-на-n матрицу нулей. (spzeros(m,n) пуста.)
spones(S) ones(m,n) Создаёт матрицу, заполненную единицами. В отличие от плотной версии, spones() имеет ту же структуру разреженности, что и S.
speye(n) eye(n) Создаёт n-на-n единичную матрицу.
full(S) sparse(A) Преобразует между плотным и разреженным форматами.
sprand(m,n,d) rand(m,n) Создаёт m-на-n случайную матрицу (плотностью d) с независимыми и одинаково распределёнными ненулевыми элементами, равномерно распределёнными на полуоткрытом интервале \([0, 1)\).
sprandn(m,n,d) randn(m,n) Создаёт m-на-n случайную матрицу (плотностью d) с независимыми и одинаково распределёнными ненулевыми элементами, распределёнными по стандартному нормальному (гауссову) распределению.
sprandn(m,n,d,X) randn(m,n,X) Создаёт m-на-n случайную матрицу (плотностью d) с независимыми и одинаково распределёнными ненулевыми элементами, распределёнными по распределению X. (Требуется пакет Distributions.)

© 2009–2016 Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors
Licensed under the MIT License.
https://docs.julialang.org/en/release-0.5/manual/arrays/

Spec-Zone.ru

Настройки Оффлайн Что нового Помощь О нас
Spec-Zone .ru
спецификации, руководства, описания, API