Juliaでは,スカラと配列は明確に区別されます。たとえばスカラと,$1 \times 1$の2次元配列は全く異なるものとして扱われます(これ以降は,1次元配列をベクトル,2次元配列を行列と表記します)。線形代数を多用する場面では,与えた式が次元に関して整合するなら,それらの区別が問題になることはありません(問題にならないようにプログラムを作るべきだ,という言い方もできます)。しかし,幸か不幸か,スカラ・ベクトル・行列を相互変換したいケースが発生するかもしれません。ここでは,それらの変換について説明します。
なお,ここでは,要素が実数であるケースのみを扱います。また、最大で二次元配列(行列)までのみを対象にします。
スカラをベクトルにする
スカラを []
で囲むと,ベクトル(列ベクトル)が得られます。
[3.0] # 列ベクトル,つまりArray{Int64,1}になる
また,関数fill()
は,特定の値(スカラ)で埋めた任意の次元の配列を作ります。
fill(3.0,1) # 要素が1個の列ベクトルになる
なお,転置記号を使えば,行ベクトルになります。
[3.0]' # 行ベクトル,つまりLinearAlgebra.Adjoint{Int64,Array{Int64,1}}になる
transpose([3.0]) # 転置ベクトル,つまりLinearAlgebra.Transpose{Int64,Array{Int64,1}}になる
[3.0]'==transpose([3.0]) # この場合の両者は等価である
ベクトルをスカラにする
ベクトルの要素を[]
で切り出します。
a = [1,2,3] # 列ベクトル
a[1] # 1番目の要素をスカラとして扱う
b = [1,2,3]' # 行ベクトル
b[1] # 1番目の要素をスカラとして扱う
スカラを行列にする
上で紹介したfill()
で1行1列の行列を作ることができます。
fill(3.0,1,1) # 1x1の行列を作る
いちどスカラをベクトルに変換し,reshape()
で1行1列の行列に変換してもOKです。
reshape([3.0],1,1) # 1x1の行列を作る
すべての要素が1の行列を作り,それに任意のスカラを乗じてもよいでしょう。
ones(1,1)*3.0 # 1x1の行列を作る
上記の関数が思い出せないときは,いちど1行1列の行列を作って変数に保存し,その要素を書き換えたほうがわかりやすいでしょう。
M = Matrix(undef,1,1) # 1x1の行列を作る
M[1,1] = 3.0
行列をスカラにする
行列の要素を[]
で切り出します。
M = [1.0 2.0; 4.0 5.0]
M[1,1] # 1,1の要素をスカラとして扱う
ベクトルを行列にする
おそらく汎用性が高いのはreshape()
関数でしょう。
a = [1,2,3] # 列ベクトル
reshape(a,3,1) # 3行1列の行列 Array{Int64,2}
reshape(a,1,3) # 1行3列の行列 Array{Int64,2}
行ベクトルをreshapeする場合,生成される行列は通常の2次元配列ではありません。このままでも演算に利用できるはずですが,少々複雑な型になっています。
b = [1,2,3]' # 行ベクトル
reshape(b,3,1) # 3行1列の行列 reshape(::LinearAlgebra.Adjoint{Int64,Array{Int64,1}}, 3, 1) with eltype Int64
reshape(b,1,3) # 1行3列の行列 reshape(::LinearAlgebra.Adjoint{Int64,Array{Int64,1}}, 1, 3) with eltype Int64
c = transpose([1,2,3])
reshape(c,3,1) # 3行1列の行列 reshape(::LinearAlgebra.Transpose{Int64,Array{Int64,1}}, 1, 3) with eltype Int64
reshape(c,1,3) # 1行3列の行列 reshape(::LinearAlgebra.Adjoint{Int64,Array{Int64,1}}, 3, 1) with eltype Int64
生成された行列をcollect
することで,通常の行列に変換できます。
b = [1,2,3]'
collect(reshape(b,3,1)) # 3行1列の行列 Array{Int64,2}
collect(reshape(b,1,3)) # 1行3列の行列 Array{Int64,2}
c = transpose([1,2,3])
collect(reshape(c,3,1)) # 3行1列の行列 Array{Int64,2}
collect(reshape(c,1,3)) # 1行3列の行列 Array{Int64,2}
あるいは,cat
,hcat
,vcat
を使えば,変換後の要素数を指定することなく,ベクトルから行列に一気に変換できます。ただし使い方を誤ると,ベクトルのままになります。
a = [1,2,3]
vcat(a) # ベクトルのまま
hcat(a) # 行列への変換
cat(a,dims=2) # 行列への変換
b = [1,2,3]'
vcat(b) # 行列への変換
hcat(b) # ベクトルのまま
cat(b,dims=1) # 行列への変換
cat(b,dims=2) # 同じく行列への変換
c = transpose([1,2,3])
vcat(c) # 行列への変換
hcat(c) # ベクトルのまま
cat(c,dims=1) # 行列への変換
cat(c,dims=2) # 同じく行列への変換
おそらく,この種の変換方法は数多くあるはずで,利用者の好みで使い分けることになるでしょう。
配列の配列とベクトル・行列の変換
配列の配列(ベクトルまたは行列のベクトル)を行列に展開する
ベクトルのベクトルは、hcat()
またはvcat()
で行列に変換できます。
julia> a=[[1,2], [3,4], [5,6]] # 2要素のベクトルの配列
julia> hcat(a...) # 水平方向に連結=行列
2×3 Array{Int64,2}:
1 3 5
2 4 6
julia> vcat(a...) # 垂直方向に連結=ベクトル
6-element Array{Int64,1}:
1
2
3
4
5
6
同様に、行列のベクトルも、ひとつの行列に展開できます。
julia> a=[[1 2; 2 1], [3 4; 4 3], [5 6; 6 5]] # 2x2行列のベクトル
3-element Array{Array{Int64,2},1}:
[1 2; 2 1]
[3 4; 4 3]
[5 6; 6 5]
julia> hcat(a...) # 水平方向に連結=行列
2×6 Array{Int64,2}:
1 2 3 4 5 6
2 1 4 3 6 5
julia> vcat(a...) # 垂直方向に連結=行列
6×2 Array{Int64,2}:
1 2
2 1
3 4
4 3
5 6
6 5
行列を配列の配列にする
内包表現を使えば可能です。上の最初の例の逆変換です。
julia> B=[
1 3 5
2 4 6
]
2×3 Array{Int64,2}:
1 3 5
2 4 6
julia> [B[:,i] for i in 1:3]
3-element Array{Array{Int64,1},1}:
[1, 2]
[3, 4]
[5, 6]
行列に拡張した例です。
julia> C=[
1 2 3 4 5 6
2 1 4 3 6 5
]
2×6 Array{Int64,2}:
1 2 3 4 5 6
2 1 4 3 6 5
julia> [C[:,i:i+1] for i in 1:2:6]
3-element Array{Array{Int64,2},1}:
[1 2; 2 1]
[3 4; 4 3]
[5 6; 6 5]
行列をベクトルにする
ここでは,行列全体をベクトルに変換する方法と,行列の一部をベクトルとして扱う方法を紹介します。例によって多くの方法があるので,そのうちの一部だけを紹介します。
行列全体をベクトルに変換する
関数vec()
が利用できます。これは,REPLでベクトルのつもりが誤って行列を作ってしまった場合の変換にも役立ちます。
M = Matrix(collect(1:12),4,3)
vec(M) # 列を縦方向につなげたベクトルを作る
# ベクトルのつもりがカンマを忘れて行列になってしまった場合の救済策
vec([1 2 3 4 5])
あるいは,行列の添字を単に[:]
とすれば,vec()
と同じ効果が得られます。
M = Matrix(collect(1:12),4,3)
M[:] # vec(M)と同じ
# vec([1 2 3 4 5])と同じ
[1 2 3 4 5][:]
ところで,特に対称行列を扱う場合,上三角または下三角の要素のみをベクトルとして抜き出したいことがあります。線形代数ではvech
という関数で呼ばれることがあります。Juliaにはこの関数は存在しませんが,これを実現するのは容易です。
A = [3 2 1; 2 4 1; 1 1 5]
A[tril(trues(size(A)))] # 下三角要素をベクトルにする
A[triu(trues(size(A)))] # 上三角要素をベクトルにする
# 関数を定義してもよい
vech(A) = A[tril(trues(size(A)))]
行列の一部をベクトル(または行列)として扱う
ブラケット[]
で,特定の行または列を抜き出します。このように抜き出した行または列は,ベクトルとして扱われます。この方法は,行列から行列を抜き出す場合にも有効です。
M = Matrix(collect(1:12),4,3)
M[:,1] # 1列目だけを抜き出す
部分配列をビューとして保持する
さて,上記のように配列の一部の要素を抜き出すとき(部分配列またはスライスを作るとき),Juliaではその全要素がコピーされ,新たなメモリ上に配置されます。そのベクトルをプログラム内で使い回すなら,テンポラリメモリへのコピーは許容できるかもしれません。しかし,単に行列の一部分を一時的に参照したいだけなら,このコピーは無駄です。
Juliaでは,コピーを作らずに行列の一部(より一般には配列の一部)を参照するために,ビューという仕組みがあります。これは,配列への参照(特定要素のインデックス)を作る仕組みで,要素そのものをコピーすることはありません。配列への参照が頻繁に発生する場合,メモリ消費量を減らして高速化するためにビューが多用されます。
以下に簡単なベンチマークをとります。これは,行列M
の各列$j$の平方和を計算し,その和を返すものです。
# 素直なプログラム
function f1(M) a=0; for j=1:3; a = a + M[:,j]'*M[:,j]; end; return a; end
# view関数を使ったプログラム
function f2(M); a=0; for j=1:3; a = a + view(M,:,j)'*view(M,:,j); end; return a; end
# viewマクロを使ったプログラム
function f3(M); a=0; for j=1:3; a = a + (@view M[:,j])'*(@view M[:,j]); end; return a; end
# 単にループを回すプログラム
function f4(M); a=0; for j=1:3; for i=1:4; a = a + M[i,j]*M[i,j]; end; end; return a; end
2020年10月9日追記:Julia 1.5にて実施しています。エスケープ忘れも修正しました。
# ベンチマーク
julia> using BenchmarkTools
julia> M = reshape(collect(1:12),4,3)
julia> @btime f1($M)
229.257 ns (6 allocations: 672 bytes)
650
julia> @btime f2($M)
36.629 ns (0 allocations: 0 bytes)
650
julia> @btime f3($M)
37.133 ns (0 allocations: 0 bytes)
650
julia> @btime f4($M)
6.305 ns (0 allocations: 0 bytes)
650
まったく同じ処理をしているのに,一番最初のプログラムが最も遅く,最も多くのメモリを消費しています。これは,ループ内でM[:,j]
が参照されるたび,それがテンポラリメモリにコピーされているためです。
真ん中の2つは,関数またはマクロを使ってM
への参照を作成しています。Julia 1.5以降では、この場合にはメモリを消費しません(ビューはスタックに確保されています)。
一番下のプログラムは,ループを使って1要素ずつ要素を参照しています。これもメモリ消費がなく,最も高速です。