LinearAlgebra

多次元配列のサポートに加えて (およびその一部として)、Julia は一般的で有用な多くの線形代数演算に対するネイティブ実装を提供します。この実装は using LinearAlgebra で読み込むことができます。tr, det, inv のような基本的な演算は全てサポートされます:

julia> A = [1 2 3; 4 1 6; 7 8 1]
3×3 Array{Int64,2}:
 1  2  3
 4  1  6
 7  8  1

julia> tr(A)
3

julia> det(A)
104.0

julia> inv(A)
3×3 Array{Float64,2}:
 -0.451923   0.211538    0.0865385
  0.365385  -0.192308    0.0576923
  0.240385   0.0576923  -0.0673077

固有値や固有ベクトルの計算といったその他の有用な演算もサポートされます:

julia> A = [-4. -17.; 2. 2.]
2×2 Array{Float64,2}:
 -4.0  -17.0
  2.0    2.0

julia> eigvals(A)
2-element Array{Complex{Float64},1}:
 -1.0 - 5.0im
 -1.0 + 5.0im

julia> eigvecs(A)
2×2 Array{Complex{Float64},2}:
  0.945905-0.0im        0.945905+0.0im
 -0.166924+0.278207im  -0.166924-0.278207im

さらに、Julia は様々な分解演算をサポートします。事前に行列を扱いやすい (実行速度やメモリ使用量に関して有利な) 形に分解することで、線形方程式の求解や行列のべき乗といった問題を高速化できます。詳細は factorize のドキュメントを参照してください。使用例を示します:

julia> A = [1.5 2 -4; 3 -1 -6; -10 2.3 4]
3×3 Array{Float64,2}:
   1.5   2.0  -4.0
   3.0  -1.0  -6.0
 -10.0   2.3   4.0

julia> factorize(A)
LU{Float64,Array{Float64,2}}
L factor:
3×3 Array{Float64,2}:
  1.0    0.0       0.0
 -0.15   1.0       0.0
 -0.3   -0.132196  1.0
U factor:
3×3 Array{Float64,2}:
 -10.0  2.3     4.0
   0.0  2.345  -3.4
   0.0  0.0    -5.24947

この例で A はハミルトン行列・対称行列・三角行列・三重対角行列・二重対角行列のいずれでもないので、この行列に行えるのは LU 分解が精一杯でしょう。しかし次の例では異なります:

julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
3×3 Array{Float64,2}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

julia> factorize(B)
BunchKaufman{Float64,Array{Float64,2}}
D factor:
3×3 Tridiagonal{Float64,Array{Float64,1}}:
 -1.64286   0.0   
  0.0      -2.8  0.0
           0.0  5.0
U factor:
3×3 UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  0.142857  -0.8
     1.0       -0.6
               1.0
permutation:
3-element Array{Int64,1}:
 1
 2
 3

この例で Julia は B が実際には対称行列であることを検出し、LU 分解よりも適切な分解を行っています。対称行列や三重対角行列といった特別な性質を持つ行列に対しては、より効率的なコードが書ける場合が多くあります。Julia では、行列がこういった性質を持つことを表す “タグ” を付けるための特別な型が提供されています:

julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
3×3 Array{Float64,2}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

julia> sB = Symmetric(B)
3×3 Symmetric{Float64,Array{Float64,2}}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

sB には (実) 対称行列であるというタグが付いているので、今後 sB に対して行われる固有値分解や行列-ベクトル積の計算といった演算は半分の要素だけを参照する効率的なものになります:

julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
3×3 Array{Float64,2}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

julia> sB = Symmetric(B)
3×3 Symmetric{Float64,Array{Float64,2}}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

julia> x = [1; 2; 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> sB\x
3-element Array{Float64,1}:
 -1.7391304347826084
 -1.1086956521739126
 -1.4565217391304346

ここで \ は線形方程式の求解演算を表します。この左除算演算子は非常に強力であり、これを使えばあらゆる種類の線形方程式を解くことができるだけの柔軟性を持ったコンパクトで読みやすいコードを簡単に記述できます。

特殊行列

線形代数では行列が特殊な対称性や構造を持ち、それが行列の分解に関連付けられることが多くあります。Julia は特殊な行列型を豊富に提供するので、特定の行列型に対して個別に開発される特殊化されたルーチンを使った高速な計算が可能です。

Julia で実装される特殊行列型を次の表にまとめます。LAPACK の最適化されたメソッドに対するフックが利用できるかどうかも示されています:

基本的な演算

行列型 + - * \ 最適化されている他のメソッド
Symmetric MV inv, sqrt, exp
Hermitian MV inv, sqrt, exp
UpperTriangular MV MV inv, det
UnitUpperTriangular MV MV inv, det
LowerTriangular MV MV inv, det
UnitLowerTriangular MV MV inv, det
UpperHessenberg M inv, det
SymTridiagonal M M MS MV eigmax, eigmin
Tridiagonal M M MS MV
Bidiagonal M M MS MV
Diagonal M M MV MV inv, det, logdet, /
UniformScaling M M MVS MVS /

凡例:

キー 説明
M (matrix) 行列-行列演算に対する最適化されたメソッドが利用可能
V (vector) 行列-ベクトル演算に対する最適化されたメソッドが利用可能
S (scalar) 行列-スカラー演算に対する最適化されたメソッドが利用可能

行列の分解

行列型 LAPACK eigen eigvals eigvecs svd svdvals
Symmetric SY ARI
Hermitian HE ARI
UpperTriangular TR A A A
UnitUpperTriangular TR A A A
LowerTriangular TR A A A
UnitLowerTriangular TR A A A
SymTridiagonal ST A ARI AV
Tridiagonal GT
Bidiagonal BD A A
Diagonal DI A

凡例:

説明
A (all) 全ての特性値/特性ベクトルを見つける最適化メソッドが利用可能 eigvals(M)
R (range) il 番目から ih 番目までの特性値を計算する最適化されたメソッドが利用可能 eigvals(M, il, ih)
I (interval) 区間 [vl, vh] に属する特性値を計算する最適化されたメソッドが利用可能 eigvals(M, vl, vh)
V (vector) 特性値 x=[x1, x2,...] に対応する特性ベクトルを計算する最適化されたメソッドが利用可能 eigvecs(M, x)

均等スケーリング作用素

均等スケーリング作用素 UniformScaling はスカラーと恒等作用素の積 λ*I を表します。恒等作用素 I は定数として定義される UniformScaling のインスタンスです。均等スケーリング作用素は汎用的なサイズを持ち、二項演算 +, -, *, \ で使うともう一方の行列のサイズに合わせられます。これは A+IA-I では A が正方である必要があることを意味します。恒等作用素 I との乗算は (スケーリング係数が 1 であることの確認を除けば) noop であり、オーバーヘッドはほとんどありません。

UniformScaling の使用例を示します:

julia> U = UniformScaling(2);

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

julia> a + U
2×2 Array{Int64,2}:
 3  2
 3  6

julia> a * U
2×2 Array{Int64,2}:
 2  4
 6  8

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

julia> b = [1 2 3; 4 5 6]
2×3 Array{Int64,2}:
 1  2  3
 4  5  6

julia> b - U
ERROR: DimensionMismatch("matrix is not square: dimensions are (2, 3)")
Stacktrace:
[...]

同一の A と異なる μ に対して (A+μI)x = b という形をした系をいくつも解く必要がある場合には、hessenberg 関数を使って A のヘッセンベルグ分解 F を求めると効率が向上する可能性があります。F が与えられると、Julia は効率的なアルゴリズムを使って (F+μ*I) \ b (元の行列を使った (A+μ*I)x \ b と等価な演算) や行列式の計算といった演算を行います。

行列の分解

行列の分解 (factorication/decomposition) は一つの行列を複数の行列の積に分解する操作であり、線形代数における中心的な問題の一つです。

Julia で実装される行列の分解を次の表にまとめます。関連するメソッドの詳細は線形代数のドキュメントの標準的な関数の節にあります。

標準的な関数

Julia が持つ線形代数関数の多くは LAPACK の関数を呼び出すことで実装されます。疎行列の分解では SuiteSparse の関数が呼ばれます。

Base.:* ── メソッド
*(A::AbstractMatrix, B::AbstractMatrix)

行列の乗算です。

julia> [1 1; 0 1] * [1 0; 1 1]
2×2 Array{Int64,2}:
 2  1
 1  1
Base.:\ ── メソッド
\(A, B)

多重アルゴリズム1を使った行列除算です。入力の行列 AB に対して、A が正方行列のとき返り値は A*X == B を満たす X です。使われるソルバは A の構造に依存します。A が上三角行列または下三角行列 (あるいは対角行列) なら A の分解は必要とならず、系は後退代入または前進代入を使って解かれます。三角でない正方行列に対しては LU 分解が使われます。

矩形行列 A に対する返り値は最小二乗最小ノルム解です。A のピボット選択付き QR 分解と R 因子に基づいた A のランク推定を使って計算されます。

A が疎行列のときは同様の多重アルゴリズムが使われます。不定値行列 A に対する LDLt 分解は数値的分解でピボット選択を行わないので、可逆行列であっても分解手続きが失敗する可能性があります。

julia> A = [1 0; 1 -2]; B = [32; -4];

julia> X = A \ B
2-element Array{Float64,1}:
 32.0
 18.0

julia> A * X == B
true
SingularException

入力行列が一つ以上のゼロ固有値を持ち可逆でないときに送出される例外です。そのような行列が絡む線形方程式系の解は計算できません。info フィールドがゼロ固有値の位置を示します。

PosDefException

入力行列が正定値でないときに送出される例外です。線形代数の関数と分解には正定値行列だけに適用できるものが存在します。info フィールドがゼロ以下の固有値の位置を示します。

ZeroPivotException <: Exception

行列の分解や線形方程式系の求解でピボットの位置 (対角要素) にゼロが現れ手続きが進められなくなったときに送出される例外です。これは行列が特異であることを意味するとは限りません。変数を並べ替えることで疑わしいゼロピボットを取り除くことができるピボット選択付きの LU 分解などに切り替えると分解が成功する可能性があります。info フィールドがゼロピボットの位置を示します。

LinearAlgebra.dot ── 関数
dot(x, y)
x ⋅ y

二つのベクトルの間のドット積を計算します。複素ベクトルに対しては、一つ目のベクトルの共役が取られます。

dot が要素間に定義される限り、dot は任意の反復可能オブジェクト (例えば任意の次元を持つ配列) に対しても動作します。

引数が同じ長さを持つという仮定の下で、dotsum(dot(vx,vy) for (vx,vy) in zip(x, y)) と意味論的に同じです。

x ⋅ ydot(x, y) と等価です ( は REPL で \cdot とタイプしてからタブ補完すると入力できます)。

julia> dot([1; 1], [2; 3])
5

julia> dot([im; im], [1; 1])
0 - 2im

julia> dot(1:5, 2:6)
70

julia> x = fill(2., (5,5));

julia> y = fill(3., (5,5));

julia> dot(x, y)
150.0
LinearAlgebra.cross ── 関数
cross(x, y)
×(x,y)

二つの三要素ベクトルに対するクロス積を計算します。

julia> a = [0;1;0]
3-element Array{Int64,1}:
 0
 1
 0

julia> b = [0;0;1]
3-element Array{Int64,1}:
 0
 0
 1

julia> cross(a,b)
3-element Array{Int64,1}:
 1
 0
 0
LinearAlgebra.factorize ── 関数
factorize(A)

行列 A の性質に基づいて、使いやすい分解を計算します。一般の行列 A が渡されると、factorizeA が対称行列や三角行列といった性質を持つかどうかの判定を行います。factorizeA の性質を確認するために A の全ての要素を確認しますが、性質が成り立たないことが分かると残りの要素を飛ばします (短絡評価を行います)。factorize の返り値は複数の系を効率的に解くのに利用可能です。例えば A=factorize(A); x=A\b; y=A\C などとしてください。A の性質と factorize(A) が行う分解の対応は次の通りです:

A の性質 分解のタイプ
正定値 コレスキー (cholesky)
密かつ対称/エルミート バンチ-カウフマン (bunchkaufman)
疎かつ対称/エルミート LDLt (ldlt)
三角 三角
対角 対角
二重対角 二重対角
三角対角 LU (lu)
実対称三重対角 LDLt (ldlt)
一般的な正方行列 LU (lu)
一般的な矩形行列 QR (qr)

例えば正定値エルミート行列に対して factorize を呼び出すと、コレスキー分解が返ります。

julia> A = Array(Bidiagonal(fill(1.0, (5, 5)), :U))
5×5 Array{Float64,2}:
 1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0
 0.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  0.0  1.0

julia> factorize(A) # factorize は A が最初から分解になっているかどうかを確認する。
5×5 Bidiagonal{Float64,Array{Float64,1}}:
 1.0  1.0           
     1.0  1.0       
         1.0  1.0   
             1.0  1.0
                 1.0

この例で factorize が返しているのは 5×5 Bidiagonal{Float64} 型の値であり、これを他の線形代数関数 (例えば固有値ソルバ) に渡すと Bidiagonal 型に対して特殊化されたメソッドが使われます。

Diagonal(A::AbstractMatrix)

A の対角要素から対角行列を構築します。

julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> Diagonal(A)
3×3 Diagonal{Int64,Array{Int64,1}}:
 1    
   5  
     9
Diagonal(V::AbstractVector)

V を対角要素とする対角行列を構築します。

julia> V = [1, 2]
2-element Array{Int64,1}:
 1
 2

julia> Diagonal(V)
2×2 Diagonal{Int64,Array{Int64,1}}:
 1  
   2
Bidiagonal(dv::V, ev::V, uplo::Symbol) where V <: AbstractVector

上 (uplo=:U) または下 (uplo=:L) の二重対角行列を対角要素 dv と副対角要素 ev のベクトルから構築します。返り値は Bidiagonal 型であり、この型に関して特殊化された効率的な線形ソルバが提供されます。通常の配列への変換は convert(Array, _) (または省略形 Array(_)) で行えます。ev の長さは dv の長さより 1 だけ小さい必要があります。

julia> dv = [1, 2, 3, 4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> ev = [7, 8, 9]
3-element Array{Int64,1}:
 7
 8
 9

julia> Bu = Bidiagonal(dv, ev, :U) # ev は対角要素の一つ上に並ぶ
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1  7    
   2  8  
     3  9
       4

julia> Bl = Bidiagonal(dv, ev, :L) # ev は対角要素の一つ下に並ぶ
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1      
 7  2    
   8  3  
     9  4
Bidiagonal(A, uplo::Symbol)

行列 A から Bidiagonal 行列を構築します。A の対角要素と、A の第一優対角要素 (uplo=:U のとき) または第一劣対角要素 (uplo=:L のとき) が使われます。

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

julia> Bidiagonal(A, :U) # A の主対角要素と第一優対角要素からなる二重対角行列
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1  1    
   2  2  
     3  3
       4

julia> Bidiagonal(A, :L) # A の主対角要素と第一劣対角要素からなる二重対角行列
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1      
 2  2    
   3  3  
     4  4
SymTridiagonal(dv::V, ev::V) where V <: AbstractVector

対称三重対角行列を対角要素 dv と第一優/劣対角要素 ev から構築します。返り値は SymTridiagonal 型であり、この型に関して特殊化された効率的な固有値ソルバが提供されます。通常の行列への変換は convert(Array, _) (または省略形 Array(_)) で行えます。

SymTridiagonal なブロック行列では dv の要素が対称化され、引数 ev は優対角ブロックとみなされます。そのとき劣対角ブロックは対応する優対角ブロックの (コピーされた) 転置となります。

julia> dv = [1, 2, 3, 4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> ev = [7, 8, 9]
3-element Array{Int64,1}:
 7
 8
 9

julia> SymTridiagonal(dv, ev)
4×4 SymTridiagonal{Int64,Array{Int64,1}}:
 1  7    
 7  2  8  
   8  3  9
     9  4

julia> A = SymTridiagonal(fill([1 2; 3 4], 3), fill([1 2; 3 4], 2));

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

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

julia> A[2,1]
2×2 Array{Int64,2}:
 1  3
 2  4
SymTridiagonal(A::AbstractMatrix)

対称行列 A の対角要素と第一優対角要素から対称三重対角行列を構築します。

julia> A = [1 2 3; 2 4 5; 3 5 6]
3×3 Array{Int64,2}:
 1  2  3
 2  4  5
 3  5  6

julia> SymTridiagonal(A)
3×3 SymTridiagonal{Int64,Array{Int64,1}}:
 1  2  
 2  4  5
   5  6

julia> B = reshape([[1 2; 2 3], [1 2; 3 4], [1 3; 2 4], [1 2; 2 3]], 2, 2);

julia> SymTridiagonal(B)
2×2 SymTridiagonal{Array{Int64,2},Array{Array{Int64,2},1}}:
 [1 2; 2 3]  [1 3; 2 4]
 [1 2; 3 4]  [1 2; 2 3]
Tridiagonal(dl::V, d::V, du::V) where V <: AbstractVector

引数に与えられる第一劣対角要素 dl, 対角要素 d, 第一優対角要素 du から三重対角行列を構築します。返り値の型は Tridiagonal であり、この型に関して特殊化された効率的な線形ソルバが提供されます。通常の行列への変換は convert(Array, _) (または省略形 Array(_)) で行えます。dldu の長さは d の長さより 1 だけ小さくなければなりません。

julia> dl = [1, 2, 3];

julia> du = [4, 5, 6];

julia> d = [7, 8, 9, 0];

julia> Tridiagonal(dl, d, du)
4×4 Tridiagonal{Int64,Array{Int64,1}}:
 7  4    
 1  8  5  
   2  9  6
     3  0
Tridiagonal(A)

行列 A の第一劣対角要素・対角要素・第一優対角要素から三重対角行列を構築します。

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

julia> Tridiagonal(A)
4×4 Tridiagonal{Int64,Array{Int64,1}}:
 1  2    
 1  2  3  
   2  3  4
     3  4
Symmetric(A, uplo=:U)

行列 A の上 (uplo = :U) または下 (uplo = :L) 三角要素に対するビューを使って Symmetric を構築します。

julia> A = [1 0 2 0 3; 0 4 0 5 0; 6 0 7 0 8; 0 9 0 1 0; 2 0 3 0 4]
5×5 Array{Int64,2}:
 1  0  2  0  3
 0  4  0  5  0
 6  0  7  0  8
 0  9  0  1  0
 2  0  3  0  4

julia> Supper = Symmetric(A)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  2  0  3
 0  4  0  5  0
 2  0  7  0  8
 0  5  0  1  0
 3  0  8  0  4

julia> Slower = Symmetric(A, :L)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  6  0  2
 0  4  0  9  0
 6  0  7  0  3
 0  9  0  1  0
 2  0  3  0  4

A 自身が対称 (つまり A == transpose(A)) でない限り SupperSlower は同じにならないことに注意してください。

Hermitian(A, uplo=:U)

行列 A の上 (uplo = :U) または下 (uplo = :L) 三角要素に対するビューを使って Hermitian を構築します。

julia> A = [ 1+0im  0+0im  2+2im  0+0im  3-3im
             0+0im  4+0im  0+0im  5+0im  0+0im
             6-6im  0+0im  7+0im  0+0im  8+8im
             0+0im  9+0im  0+0im  1+0im  0+0im
             2+2im  0+0im  3-3im  0+0im  4+0im ];

julia> Hupper = Hermitian(A)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  2+2im  0+0im  3-3im
 0+0im  4+0im  0+0im  5+0im  0+0im
 2-2im  0+0im  7+0im  0+0im  8+8im
 0+0im  5+0im  0+0im  1+0im  0+0im
 3+3im  0+0im  8-8im  0+0im  4+0im

julia> Hlower = Hermitian(A, :L)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  6+6im  0+0im  2-2im
 0+0im  4+0im  0+0im  9+0im  0+0im
 6-6im  0+0im  7+0im  0+0im  3+3im
 0+0im  9+0im  0+0im  1+0im  0+0im
 2+2im  0+0im  3-3im  0+0im  4+0im

A 自身がエルミート (つまり A == adjoint(A)) でない限り HupperHlower は同じにならないことに注意してください。

対角成分の虚部は無視されます。

Hermitian(fill(complex(1,1), 1, 1)) == fill(1, 1, 1)
LowerTriangular(A::AbstractMatrix)

行列 A に対するビューを使って LowerTriangular を構築します。

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  9.0

julia> LowerTriangular(A)
3×3 LowerTriangular{Float64,Array{Float64,2}}:
 1.0       
 4.0  5.0   
 7.0  8.0  9.0
UpperTriangular(A::AbstractMatrix)

行列 A に対するビューを使って UpperTriangular を構築します。

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  9.0

julia> UpperTriangular(A)
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 1.0  2.0  3.0
     5.0  6.0
         9.0
UnitLowerTriangular(A::AbstractMatrix)

行列 A に対するビューを使って UnitLowerTriangular を構築します。返り値は対角要素に Aeltypeoneunit を持ちます。

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  9.0

julia> UnitLowerTriangular(A)
3×3 UnitLowerTriangular{Float64,Array{Float64,2}}:
 1.0       
 4.0  1.0   
 7.0  8.0  1.0
UnitUpperTriangular(A::AbstractMatrix)

行列 A に対するビューを使って UnitUpperTriangular を構築します。返り値は対角要素に Aeltypeoneunit を持ちます。

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  9.0

julia> UnitUpperTriangular(A)
3×3 UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  2.0  3.0
     1.0  6.0
         1.0
UpperHessenberg(A::AbstractMatrix)

行列 A に対するビューを使って UpperHessenberg を構築します。A の第一劣対角要素より下の部分は無視されます。

H \ bdet(H) に対する効率的なアルゴリズムが実装されます。

任意の行列を同様の上ヘッセンベルグ行列に分解するには hessenberg 関数を参照してください2

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

julia> UpperHessenberg(A)
4×4 UpperHessenberg{Int64,Array{Int64,2}}:
 1   2   3   4
 5   6   7   8
   10  11  12
      15  16
UniformScaling{T<:Number}

一般的なサイズを持つ均等スケーリング作用素です。均等スケーリング作用素はスカラーと恒等作用素の積 λ*I として定義されます。I も参照してください。

julia> J = UniformScaling(2.)
UniformScaling{Float64}
2.0*I

julia> A = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> J*A
2×2 Array{Float64,2}:
 2.0  4.0
 6.0  8.0
LinearAlgebra.I ── 定数
I

UniformScaling 型のオブジェクトです。任意のサイズの恒等行列を表します。

julia> fill(1, (5,6)) * I == fill(1, (5,6))
true

julia> [1 2im 3; 1im 2 3] * I
2×3 Array{Complex{Int64},2}:
 1+0im  0+2im  3+0im
 0+1im  2+0im  3+0im
LinearAlgebra.Factorization

行列の分解を表す抽象型です。利用可能な行列の分解のリストについてはドキュメントを参照してください。

LinearAlgebra.LU ──
LU <: Factorization

正方行列の LU 分解を表す型です。対応する行列分解関数 lu はこの型の値を返します。

分解 F::LU の個別要素に対するアクセスには getproperty を使います:

要素 説明
F.L LU の L (単位下三角行列)
F.U LU の U (上三角行列)
F.p (右) 置換ベクトル
F.P (右) 置換行列

LU 型の分解オブジェクトを反復すると F.L, F.U, F.p が生成されます。

julia> A = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3

julia> F = lu(A)
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.666667  1.0
U factor:
2×2 Array{Float64,2}:
 6.0  3.0
 0.0  1.0

julia> F.L * F.U == A[F.p, :]
true

julia> l, u, p = lu(A); # 反復を使った分解の展開

julia> l == F.L && u == F.U && p == F.p
true
LinearAlgebra.lu ── 関数
lu(A, pivot=Val(true); check = true) -> F::LU

A の LU 分解を計算します。

check = true だと分解が失敗したときに例外が送出されます。check = false の場合は、分解の正当性を確認する (issuccess を呼び出す) 責任はユーザーにあります。

多くの場合で、AAbstractMatrix{T} の部分型 ST 型が +, -, *, / をサポートするなら、返り値の型は LU{T,S{T}} となります。ピボット選択が有効なときは、要素型は abs< をサポートするべきです。

分解 F の個別の要素には getproperty でアクセスできます:

要素 説明
F.L LU 分解の L (単位下三角行列)
F.U LU 分解の U (上三角行列)
F.p (右) 置換ベクトル
F.P (右) 置換行列

LU 型の分解オブジェクトを反復すると F.L, F.U, F.p が生成されます。

FA の関係は次の通りです:

F.L*F.U == A[F.p, :]

F はさらに次の関数をサポートします:

サポートされる関数 LU LU{T,Tridiagonal{T}}
/
\
inv
det
logdet
logabsdet
size

julia> A = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3

julia> F = lu(A)
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.666667  1.0
U factor:
2×2 Array{Float64,2}:
 6.0  3.0
 0.0  1.0

julia> F.L * F.U == A[F.p, :]
true

julia> l, u, p = lu(A); # 反復を使った分解の展開

julia> l == F.L && u == F.U && p == F.p
true
LinearAlgebra.lu! ── 関数
lu!(A, pivot=Val(true); check = true) -> LU

lu! は基本的に lu と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。A の要素型 (例えば整数型) で表現できない数値が分解で生じると InexactError 例外が発生します。

julia> A = [4. 3.; 6. 3.]
2×2 Array{Float64,2}:
 4.0  3.0
 6.0  3.0

julia> F = lu!(A)
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.666667  1.0
U factor:
2×2 Array{Float64,2}:
 6.0  3.0
 0.0  1.0

julia> iA = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3

julia> lu!(iA)
ERROR: InexactError: Int64(0.6666666666666666)
Stacktrace:
[...]
Cholesky <: Factorization

密な対称/エルミート正定値行列 A のコレスキー分解を表す型です。対応する行列分解関数 cholesky はこの型の値を返します。

分解 F::Cholesky の三角コレスキー因子は F.L および F.U で取得できます。

julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
3×3 Array{Float64,2}:
   4.0   12.0  -16.0
  12.0   37.0  -43.0
 -16.0  -43.0   98.0

julia> C = cholesky(A)
Cholesky{Float64,Array{Float64,2}}
U factor:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
     1.0   5.0
          3.0

julia> C.U
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
     1.0   5.0
          3.0

julia> C.L
3×3 LowerTriangular{Float64,Array{Float64,2}}:
  2.0       
  6.0  1.0   
 -8.0  5.0  3.0

julia> C.L * C.U == A
true
CholeskyPivoted

密な対称/エルミート半正定値行列 A のピボット選択付きコレスキー分解を表す型です。対応する行列分解関数 cholesky(A, Val(true)) はこの型の値を返します。

分解 F::CholeskyPivoted の三角コレスキー因子は F.L および F.U で取得できます。

julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
3×3 Array{Float64,2}:
   4.0   12.0  -16.0
  12.0   37.0  -43.0
 -16.0  -43.0   98.0

julia> C = cholesky(A, Val(true))
CholeskyPivoted{Float64,Array{Float64,2}}
U factor with rank 3:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 9.89949  -4.34366  -1.61624
          4.25825   1.1694
                   0.142334
permutation:
3-element Array{Int64,1}:
 3
 2
 1
LinearAlgebra.cholesky ── 関数
cholesky(A, Val(false); check = true) -> Cholesky

対称/エルミートの正定値密行列 A のコレスキー分解を計算し、分解を表す Cholesky 型の値を返します。行列 ASymmetric または HermitianStridedMatrix もしくは完全に対称またはエルミートな StridedMatrix である必要があります。分解 F の三角コレスキー因子は F.LF.U で取得できます。Cholesky オブジェクトに対しては関数 size, \, inv, det, logdet, isposdef が利用できます。

構築で生じた丸め誤差によりぎりぎりエルミートでない行列 A をコレスキー分解するときは、Acholesky に渡す前に Hermitian(A) とラップすることで完全にエルミートな行列として扱われるようにしてください。

check = true だと分解が失敗したときに例外が送出されます。check = false の場合は、分解の正当性を確認する (issuccess を呼び出す) 責任はユーザーにあります。

julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
3×3 Array{Float64,2}:
   4.0   12.0  -16.0
  12.0   37.0  -43.0
 -16.0  -43.0   98.0

julia> C = cholesky(A)
Cholesky{Float64,Array{Float64,2}}
U factor:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
     1.0   5.0
          3.0

julia> C.U
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
     1.0   5.0
          3.0

julia> C.L
3×3 LowerTriangular{Float64,Array{Float64,2}}:
  2.0       
  6.0  1.0   
 -8.0  5.0  3.0

julia> C.L * C.U == A
true
cholesky(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted

対称/エルミートの半正定値密行列 A のピボット選択付きコレスキー分解を計算し、分解を表す CholeskyPivoted 型の値を返します。行列 ASymmetric または HermitianStridedMatrix もしくは完全に対称またはエルミートな StridedMatrix である必要があります。分解 F の三角コレスキー因子は F.LF.U で取得できます。Cholesky オブジェクトに対しては関数 size, \, inv, det, rank が利用できます。 引数 tol はランクを決定するときの許容範囲を定めます。負の値を設定すると、マシンの精度が許容範囲となります。

構築で生じた丸め誤差によりぎりぎりエルミートでない行列 A をコレスキー分解するときは、Acholesky に渡す前に Hermitian(A) とラップすることで完全にエルミートな行列として扱われるようにしてください。

check = true だと分解が失敗したときに例外が送出されます。check = false の場合は、分解の正当性を確認する (issuccess を呼び出す) 責任はユーザーにあります。

LinearAlgebra.cholesky! ── 関数
cholesky!(A, Val(false); check = true) -> Cholesky

基本的に cholesky と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。A の要素型 (例えば整数型) で表現できない数値が分解で生じると InexactError 例外が発生します。

julia> A = [1 2; 2 50]
2×2 Array{Int64,2}:
 1   2
 2  50

julia> cholesky!(A)
ERROR: InexactError: Int64(6.782329983125268)
Stacktrace:
[...]
cholesky!(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted

基本的に cholesky と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。A の要素型 (例えば整数型) で表現できない数値が分解で生じると InexactError 例外が発生します。

lowrankupdate(C::Cholesky, v::StridedVector) -> CC::Cholesky

コレスキー分解 C をベクトル v で更新します。A = C.U'C.U のとき CC = cholesky(C.U'C.U + v*v') が成り立ちますが、CC の計算で行われる演算は O(n^2) 回だけで済みます。

lowrankdowndate(C::Cholesky, v::StridedVector) -> CC::Cholesky

コレスキー分解 C をベクトル v で逆更新 (downdate) します。A = C.U'C.U のとき CC = cholesky(C.U'C.U - v*v') が成り立ちますが、CC の計算で行われる演算は O(n^2) 回だけで済みます。

lowrankupdate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

コレスキー分解 C をベクトル v で更新します。A = C.U'C.U のとき CC = cholesky(C.U'C.U + v*v') が成り立ちますが、CC の計算で行われる演算は O(n^2) 回だけで済みます。入力の分解 C はインプレースに更新されるので、関数が返ったとき C == CC が成り立ちます。ベクトル v は計算中に破壊されます。

lowrankdowndate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

コレスキー分解 C をベクトル v で逆更新 (downdate) します。A = C.U'C.U のとき CC = cholesky(C.U'C.U - v*v') が成り立ちますが、CC の計算で行われる演算は O(n^2) 回だけで済みます。入力の分解 C はインプレースに更新されるので、関数が返ったとき C == CC が成り立ちます。ベクトル v は計算中に破壊されます。

LinearAlgebra.LDLt ──
LDLt <: Factorization

SymTridiagonal 型の実対称三重対角行列 S の LDLt 分解を表す型です。LDLt 分解は S = L*Diagonal(d)*L' を満たす UnitLowerTriangular 行列 L とベクトル d からなります。LDLt 分解 F = ldlt(S) は主に線形方程式系 Sx = bF\b として解くのに使われます。対応する行列分解関数 ldlt はこの型の値を返します。

分解 F::LDLt の個別要素には getproperty を通してアクセスできます:

要素 説明
F.L LDLt 分解の L (単位下三角行列)
F.D LDLt 分解の D (対角行列)
F.Lt LDLt 分解の Lt (単位上三角行列)
F.d D の対角要素 (Vector)

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 3.0  1.0   
 1.0  4.0  2.0
     2.0  5.0

julia> F = ldlt(S)
LDLt{Float64,SymTridiagonal{Float64,Array{Float64,1}}}
L factor:
3×3 UnitLowerTriangular{Float64,SymTridiagonal{Float64,Array{Float64,1}}}:
 1.0                 
 0.333333  1.0        
 0.0       0.545455  1.0
D factor:
3×3 Diagonal{Float64,Array{Float64,1}}:
 3.0           
     3.66667   
             3.90909
LinearAlgebra.ldlt ── 関数
ldlt(S::SymTridiagonal) -> LDLt

実対称三重対角行列 S の LDLt 分解を計算します。分解は S = L*Diagonal(d)*L' を満たす UnitLowerTriangular 行列 L とベクトル d からなります。LDLt 分解 F = ldlt(S) は主に線形方程式系 Sx = bF\b として解くのに使われます。

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 3.0  1.0   
 1.0  4.0  2.0
     2.0  5.0

julia> ldltS = ldlt(S);

julia> b = [6., 7., 8.];

julia> ldltS \ b
3-element Array{Float64,1}:
 1.7906976744186047
 0.627906976744186
 1.3488372093023255

julia> S \ b
3-element Array{Float64,1}:
 1.7906976744186047
 0.627906976744186
 1.3488372093023255
LinearAlgebra.ldlt! ── 関数
ldlt!(S::SymTridiagonal) -> LDLt

基本的に ldlt と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 3.0  1.0   
 1.0  4.0  2.0
     2.0  5.0

julia> ldltS = ldlt!(S);

julia> ldltS === S
false

julia> S
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 3.0       0.333333   
 0.333333  3.66667   0.545455
          0.545455  3.90909
LinearAlgebra.QR ──
QR <: Factorization

パックされた形式で格納された行列の QR 分解を表す型です。通常は qr で作成されます。\(A\)\(m×n\) 行列のとき、QR 分解は

\[ A = Q R \]

が成り立つような直交/ユニタリ行列 \(Q\) と上三角行列 \(R\) を計算します。行列 \(Q\) はハウスホルダー鏡映子 \(v_i\) と係数 \(\tau_i\) の列として格納され、次の等式が成り立ちます:

\[ Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T) \]

QR 型の分解オブジェクトを反復すると QR が生成されます。

QR 型のオブジェクトは二つのフィールドを持ちます:

  • factorsm×n 行列です:

    • 上三角部分には \(R\) の要素が含まれ、QR オブジェクト F に対して R = triu(F.factors) が成り立ちます。

    • 対角要素より下の部分には鏡映子 \(v_i\) がパックされた形式で格納され、行列 V = I + tril(F.factors, -1) の第 \(i\) 列が \(v_i\) となります。

  • τ は長さ min(m,n) のベクトルであり、係数 \(\tau_i\) が格納されます。

QRCompactWY <: Factorization

コンパクトなブロック形式で格納した行列の QR 分解を表す型です。通常 qr で作成されます。\(A\)\(m×n\) 行列のとき、QR 分解は

\[ A = Q R \]

が成り立つような直交/ユニタリ行列 \(Q\) と上三角行列 \(R\) を計算します。格納形式は QR と似ていますが、直交/ユニタリ行列 \(Q\)コンパクト WY 形式 [Schreiber1989] で格納されます。ブロックサイズが \(n_b\) のとき分解は \(m×n\) の下台形行列 \(V\) と行列 \(T = (T_1 \; T_2 \; ... \; T_{b-1} \; T_b')\) で表されます。ここで \(b = \lceil \min(m,n) / n_b \rceil\) であり、\(T\)\(j = 1, ..., b-1\) に対するサイズ \(n_b × n_b\) の上三角行列 \(T_j\)\(j = b\) に対するサイズ \(n_b × \min(m,n) - (b-1)n_b\) の上台形行列 \(T_j'\) からなります。\(T_b'\) の上側の正方部分を \(T_b\) としたとき、次の等式が成り立ちます:

\[ Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T) = \prod_{j=1}^{b} (I - V_j T_j V_j^T) \]

ここで \(V\)\(i\) 番目の列が \(v_i\) であり、[diag(T_1); diag(T_2); …; diag(T_b)]\(i\) 番目の列が \(\tau_i\) であり、\(V\) の左側から \(m×n_b\) のブロックを切り出したものが \((V_1 \; V_2 \; ... \; V_b)\) です。qrQRCompactWY 型のオブジェクトが構築されるときは、ブロックサイズは \(n_b = \min(m, n, 36)\) と定まります。

QRCompactWY 型の分解オブジェクトを反復すると QR が生成されます。

QRCompactWY 型のオブジェクトは二つのフィールドを持ちます:

  • factorsQR と同様に m×n 行列です:

    • 上三角部分には \(R\) の要素が含まれ、QRCompactWY オブジェクト F に対して等式 R = triu(F.factors) が成り立ちます。

    • 対角要素より下の部分には鏡映子 \(v_i\) がパックされた形式で格納され、等式 V = I + tril(F.factors, -1) が成り立ちます。

  • T は上述の \(n_b × \min(m,n)\) 行列です。三角行列 \(T_j\) の対角要素より下の部分は全て無視されます。

情報

この形式を古い WY 表現 [Bischof1987] と混同しないよう注意してください。

QRPivoted <: Factorization

行列に対する列のピボット選択付き QR 分解をパックされた形式で格納したものを表す型です。通常は qr で作成されます。\(A\)\(m×n\) 行列のとき、列のピボット選択付き QR 分解は

\[ A P = Q R \]

が成り立つような置換行列 \(P\) と直交/ユニタリ行列 \(Q\) と上三角行列 \(R\) を計算します。行列 \(Q\) はハウスホルダー鏡映子の列として格納されます:

\[ Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T) \]

QRPivoted 型の分解オブジェクトを反復すると Q, R, p が生成されます。

QRPivoted 型のオブジェクトは次のフィールドを持ちます:

  • factorsm×n 行列です:

    • 上三角部分には \(R\) の要素が格納され、QRPivoted オブジェクト F に対して R = triu(F.factors) が成り立ちます。

    • 対角要素より下の部分には鏡映子 \(v_i\) がパックされた形式で格納され、行列 V = I + tril(F.factors, -1) の第 \(i\) 列が \(v_i\) となります。

  • τ は長さ min(m,n) のベクトルで、係数 \(\tau_i\) が格納されます。

  • jpvt は長さ n の整数ベクトルで、置換 \(P\) に対応します。

LinearAlgebra.qr ── 関数
qr(A, pivot=Val(false); blocksize) -> F

行列 A の QR 分解を計算します。直交 (A が複素数の値を持つならユニタリ) 行列 Q と上三角行列 R であって次の等式を満たすものが計算されます:

\[ A = Q R \]

返り値のオブジェクト F は分解をパックされた形式で格納します:

  • pivot == Val(true) だと、FQRPivoted オブジェクトになります。

  • それ以外の場合で A の要素型が BLAS に対応した型 (Float32, Float64, ComplexF32, ComplexF64) なら、FQRCompactWY オブジェクトになります。

  • このどちらでもない場合、FQR オブジェクトになります。

分解 F の個別の要素はプロパティアクセッサで取得できます:

  • F.Q: 直交/ユニタリ行列 Q
  • F.R: 上三角行列 R
  • F.p: ピボットの置換ベクトル (QRPivoted のみ)
  • F.P: ピボットの置換行列 (QRPivoted のみ)

分解オブジェクトを反復すると Q, R, p が生成されます (p は存在する場合のみ)。

QR オブジェクトに対して利用可能な関数は inv, size, \ です。A が矩形行列だと \ は最小二乗解を計算し、解が唯一でない場合は最小ノルム解が返ります。A がフルランクでないときは、(列の) ピボット選択を行う分解は最小ノルム解を求めることが要求されます。

Q は乗算でフル/正方であるようにもフル/正方でないようにも使うことができます。言い換えると、F.Q*F.RF.Q*A は両方サポートされます。Q を通常の行列に変換するには Matrix を使ってください。この操作は “薄い” Q 因子を返します。つまり Am×nm>=n なら、Matrix(F.Q) は互いに直交する単位ベクトルを持った m×n 行列を返します。“完全な” Q 因子 (m×m 行列) を取得するには F.Q*Matrix(I,m,m) としてください。m <= n だと F.Qm×m の直交行列を生成します。

QR 分解のブロックサイズは pivot == Val(false) かつ A isa StridedMatrix{<:BlasFloat} のときキーワード引数 blocksize :: Integer で指定できます。blocksize > minimum(size(A)) のとき blocksize は無視されます。詳細は QRCompactWY を参照してください。

Julia 1.4

キーワード引数 blocksize は Julia 1.4 以降でサポートされます。

julia> A = [3.0 -6.0; 4.0 -8.0; 0.0 1.0]
3×2 Array{Float64,2}:
 3.0  -6.0
 4.0  -8.0
 0.0   1.0

julia> F = qr(A)
LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.6   0.0   0.8
 -0.8   0.0  -0.6
  0.0  -1.0   0.0
R factor:
2×2 Array{Float64,2}:
 -5.0  10.0
  0.0  -1.0

julia> F.Q * F.R == A
true
情報

qr の返り値の型が一つに定まらないのは、LAPACK が行列 QR を二つの密行列ではなくコンパクトな表現を使って格納し、ハウスホルダー基本鏡映子の積が占めるメモリ容量を最小化するためです。

LinearAlgebra.qr! ── 関数
qr!(A, pivot=Val(false); blocksize)

AStridedMatrix の部分型のとき qr! は基本的に qr と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。A の要素型 (例えば整数型) で表現できない数値が分解で生じると InexactError 例外が発生します。

Julia 1.4

キーワード引数 blocksize は Julia 1.4 以降でサポートされます。

julia> a = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> qr!(a)
LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
2×2 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.316228  -0.948683
 -0.948683   0.316228
R factor:
2×2 Array{Float64,2}:
 -3.16228  -4.42719
  0.0      -0.632456

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

julia> qr!(a)
ERROR: InexactError: Int64(-3.1622776601683795)
Stacktrace:
[...]
LinearAlgebra.LQ ──
LQ <: Factorization

行列の LQ 分解を表す型です。行列 A の LQ 分解は transpose(A)QR 分解です。対応する行列分解関数 lq はこの型の値を返します。

分解オブジェクト S::LQ の下三角要素は S.L で、対角/ユニタリ要素は S.Q で取得でき、A ≈ S.L*S.Q が成り立ちます。

LQ 型の分解オブジェクトを反復すると S.LS.Q が生成されます。

julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> S = lq(A)
LQ{Float64,Array{Float64,2}} with factors L and Q:
[-8.60233 0.0; 4.41741 -0.697486]
[-0.581238 -0.813733; -0.813733 0.581238]

julia> S.L * S.Q
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> l, q = S; # 反復を使った分解の展開

julia> l == S.L &&  q == S.Q
true
LinearAlgebra.lq ── 関数
lq(A) -> S::LQ

A の LQ 分解を計算します。LQ 型の LQ 分解オブジェクト S の下三角要素は S.L で、直交/ユニタリ要素は S.Q で取得でき、A ≈ S.L*S.Q が成り立ちます。

LQ 型の分解オブジェクトを反復すると S.LS.Q が生成されます。

A の LQ 分解は transpose(A) の QR 分解であり、劣決定線形方程式 (A の列が行より多く、かつ A が行フルランクである場合) の最小ノルム解を (lq(A) \ b で) 計算するのに有用です。

julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> S = lq(A)
LQ{Float64,Array{Float64,2}} with factors L and Q:
[-8.60233 0.0; 4.41741 -0.697486]
[-0.581238 -0.813733; -0.813733 0.581238]

julia> S.L * S.Q
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> l, q = S; # 反復を使った分解の展開

julia> l == S.L &&  q == S.Q
true
LinearAlgebra.lq! ── 関数
lq!(A) -> LQ

A を作業場所として使いながら、ALQ 分解を計算します。lq も参照してください。

BunchKaufman <: Factorization

対称/エルミート行列 A のバンチ-カウフマン分解 P'UDU'P または P'LDL'P を表す型です。どちらの形が使われるかは A に上三角要素と下三角要素のどちらが保存されているかに応じて決まります (デフォルトは上三角です)。A が複素対称行列なら U'L' は共役を取らない転置 (つまり transpose(U)transpose(L)) を表します。対応する行列分解関数 bunchkaufman はこの型の値を返します。

S::BunchKaufman が分解オブジェクトなら、分解の要素は S.D, S.uplo の値に応じて S.U または S.L, S.p で取得できます。

分解オブジェクトを反復すると S.D, S.uplo の値に応じて S.U または S.L, S.p が生成されます。

julia> A = [1 2; 2 3]
2×2 Array{Int64,2}:
 1  2
 2  3

julia> S = bunchkaufman(A) # A は内部で Symmetric(A) にラップされる。
BunchKaufman{Float64,Array{Float64,2}}
D factor:
2×2 Tridiagonal{Float64,Array{Float64,1}}:
 -0.333333  0.0
  0.0       3.0
U factor:
2×2 UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  0.666667
     1.0
permutation:
2-element Array{Int64,1}:
 1
 2

julia> d, u, p = S; # 反復を使った分解の展開

julia> d == S.D && u == S.U && p == S.p
true

julia> S = bunchkaufman(Symmetric(A, :L))
BunchKaufman{Float64,Array{Float64,2}}
D factor:
2×2 Tridiagonal{Float64,Array{Float64,1}}:
 3.0   0.0
 0.0  -0.333333
L factor:
2×2 UnitLowerTriangular{Float64,Array{Float64,2}}:
 1.0        
 0.666667  1.0
permutation:
2-element Array{Int64,1}:
 2
 1
bunchkaufman(A, rook::Bool=false; check = true) -> S::BunchKaufman

対称/エルミート行列 A のバンチ-カウフマン分解 [Bunch1977] P'*U*D*U'*P または P'*L*D*L'*P を計算します。どちらの形が使われるかは A に上三角要素と下三角要素のどちらが保存されているかに応じて決まります (デフォルトは上三角です)。返り値は BunchKaufman 型のオブジェクトです。A が複素対称行列なら U'L' は共役を取らない転置 (つまり transpose(U)transpose(L)) を表すことに注意してください。

分解オブジェクトを反復すると S.D, S.uplo の値に応じて S.U または S.L, S.p が生成されます。

rooktrue だと rook pivoting が使われます。truefalse のときは使われません。

checktrue だと分解が失敗したときに例外が送出されます。checkfalse の場合は、分解の正当性を確認する (issuccess を呼び出す) 責任はユーザーにあります。

BunchKaufman オブジェクトに対して利用可能な関数は size, \, inv, issymmetric, ishermitian, getindex です。

julia> A = [1 2; 2 3]
2×2 Array{Int64,2}:
 1  2
 2  3

julia> S = bunchkaufman(A) # A は内部で Symmetric(A) にラップされる。
BunchKaufman{Float64,Array{Float64,2}}
D factor:
2×2 Tridiagonal{Float64,Array{Float64,1}}:
 -0.333333  0.0
  0.0       3.0
U factor:
2×2 UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  0.666667
     1.0
permutation:
2-element Array{Int64,1}:
 1
 2

julia> d, u, p = S; # 反復を使った分解の展開

julia> d == S.D && u == S.U && p == S.p
true

julia> S = bunchkaufman(Symmetric(A, :L))
BunchKaufman{Float64,Array{Float64,2}}
D factor:
2×2 Tridiagonal{Float64,Array{Float64,1}}:
 3.0   0.0
 0.0  -0.333333
L factor:
2×2 UnitLowerTriangular{Float64,Array{Float64,2}}:
 1.0        
 0.666667  1.0
permutation:
2-element Array{Int64,1}:
 2
 1
bunchkaufman!(A, rook::Bool=false; check = true) -> BunchKaufman

bunchkaufman! は基本的に bunchkaufman と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。

LinearAlgebra.Eigen ──
Eigen <: Factorization

正方行列の固有値分解/スペクトル分解を表す型です。対応する行列分解関数 eigen はこの型の値を返します。

F::Eigen が分解オブジェクトなら、固有値は F.values で、固有ベクトルは行列 F.vectors の列として取得できます (k 番目の固有ベクトルはスライス F.vectors[:, k] です)。

分解オブジェクトを反復すると F.valuesF.vectors が生成されます。

julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
values:
3-element Array{Float64,1}:
  1.0
  3.0
 18.0
vectors:
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> F.values
3-element Array{Float64,1}:
  1.0
  3.0
 18.0

julia> F.vectors
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> vals, vecs = F; # 反復を使った分解の展開

julia> vals == F.values && vecs == F.vectors
true
GeneralizedEigen <: Factorization

二つの行列 A, B の一般化固有値分解/一般化スペクトル分解を表す型です。対応する行列分解関数 eigen に二つの行列を引数として渡すと、この型の値が返ります。

F::GeneralizedEigen が分解オブジェクトなら、固有値は F.values で、固有ベクトルは行列 F.vectors の列として取得できます (k 番目の固有ベクトルはスライス F.vectors[:, k] です)。

GeneralizedEigen 型の分解オブジェクトを反復すると F.valuesF.vectors が生成されます。

julia> A = [1 0; 0 -1]
2×2 Array{Int64,2}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Array{Int64,2}:
 0  1
 1  0

julia> F = eigen(A, B)
GeneralizedEigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}
values:
2-element Array{Complex{Float64},1}:
 0.0 - 1.0im
 0.0 + 1.0im
vectors:
2×2 Array{Complex{Float64},2}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im

julia> F.values
2-element Array{Complex{Float64},1}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> F.vectors
2×2 Array{Complex{Float64},2}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im

julia> vals, vecs = F; # 反復を使った分解の展開

julia> vals == F.values && vecs == F.vectors
true
LinearAlgebra.eigvals ── 関数
eigvals(A; permute::Bool=true, scale::Bool=true, sortby) -> values

A の固有値を返します。

対称でない一般の行列に対しては、固有値の計算を始める前に行列をどの程度調整するかを指定できます。キーワード引数 permute, scale, sortbyeigen! と同様です。

julia> diag_matrix = [1 0; 0 4]
2×2 Array{Int64,2}:
 1  0
 0  4

julia> eigvals(diag_matrix)
2-element Array{Float64,1}:
 1.0
 4.0
eigvals(x::Number; kwargs...)

eigvals はスカラーの入力に対してスカラーを返します。

julia> eigvals(-2)
-2
eigvals(A, B) -> values

AB の一般化固有値を計算します。

julia> A = [1 0; 0 -1]
2×2 Array{Int64,2}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Array{Int64,2}:
 0  1
 1  0

julia> eigvals(A,B)
2-element Array{Complex{Float64},1}:
 0.0 - 1.0im
 0.0 + 1.0im
eigvals(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> values

A の固有値を返します。UnitRange 型の値を irange に渡せば、固有値をソートしたときの位置 (二番目から八番目など) で計算する固有値を指定できます。

julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 1.0  2.0   
 2.0  2.0  3.0
     3.0  1.0

julia> eigvals(A, 2:2)
1-element Array{Float64,1}:
 0.9999999999999996

julia> eigvals(A)
3-element Array{Float64,1}:
 -2.1400549446402604
  1.0000000000000002
  5.140054944640259
eigvals(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> values

A の固有値を返します。vlvu を渡すと、その間にある固有値だけを計算できます。

julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 1.0  2.0   
 2.0  2.0  3.0
     3.0  1.0

julia> eigvals(A, -1, 2)
1-element Array{Float64,1}:
 1.0000000000000009

julia> eigvals(A)
3-element Array{Float64,1}:
 -2.1400549446402604
  1.0000000000000002
  5.140054944640259
LinearAlgebra.eigvals! ── 関数
eigvals!(A; permute::Bool=true, scale::Bool=true, sortby) -> values

基本的に eigvals と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。キーワード引数 permute, scale, sortbyeigen と同様です。

情報

行列 A を入力として eigvals! を呼び出しても、A に固有値が格納されるわけではありません。A は作業場所として使われます。

julia> A = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> eigvals!(A)
2-element Array{Float64,1}:
 -0.3722813232690143
  5.372281323269014

julia> A
2×2 Array{Float64,2}:
 -0.372281  -1.0
  0.0        5.37228
eigvals!(A, B; sortby) -> values

基本的に eigvals と同じですが、入力 A, B のコピーを作らずに上書きすることで空間を節約します。

情報

行列 A, B を入力として eigvals! を呼び出しても、二つの行列に固有値が格納されるわけではありません。行列は作業場所として使われます。

julia> A = [1. 0.; 0. -1.]
2×2 Array{Float64,2}:
 1.0   0.0
 0.0  -1.0

julia> B = [0. 1.; 1. 0.]
2×2 Array{Float64,2}:
 0.0  1.0
 1.0  0.0

julia> eigvals!(A, B)
2-element Array{Complex{Float64},1}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> A
2×2 Array{Float64,2}:
 -0.0  -1.0
  1.0  -0.0

julia> B
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
eigvals!(A::Union{SymTridiagonal,Hermitian,Symmetric}, irange::UnitRange) -> values

基本的に eigvals と同様ですが、入力 A のコピーを作らず上書きすることで空間を節約します。irange は計算する固有値の添え字の範囲です ──例えば二番目から八番目の固有値を計算できます。

eigvals!(A::Union{SymTridiagonal,Hermitian,Symmetric}, vl::Real, vu::Real) -> values

基本的に eigvals と同様ですが、入力 A のコピーを作らず上書きすることで空間を節約します。vl は固有値を探索する範囲の下界で、vu は上界です。

LinearAlgebra.eigmax ── 関数
eigmax(A; permute::Bool=true, scale::Bool=true)

A の最大固有値を返します。オプション permute=true を指定すると上三角行列に近づけるために行列が置換され、scale=true を指定すると行列が対角要素でスケールされ行と列の絶対値が大まかに揃えられます。複素数はソートできないので、A が複素固有値を持つとき eigmax は失敗します。

julia> A = [0 im; -im 0]
2×2 Array{Complex{Int64},2}:
 0+0im  0+1im
 0-1im  0+0im

julia> eigmax(A)
1.0

julia> A = [0 im; -1 0]
2×2 Array{Complex{Int64},2}:
  0+0im  0+1im
 -1+0im  0+0im

julia> eigmax(A)
ERROR: DomainError with Complex{Int64}[0+0im 0+1im; -1+0im 0+0im]:
`A` cannot have complex eigenvalues.
Stacktrace:
[...]
LinearAlgebra.eigmin ── 関数
eigmin(A; permute::Bool=true, scale::Bool=true)

A の最小固有値を返します。オプション permute=true を指定すると上三角行列に近づけるために行列が置換され、scale=true を指定すると行列が対角要素でスケールされ行と列の絶対値が大まかに揃えられます。複素数はソートできないので、A が複素固有値を持つとき eigmin は失敗します。

julia> A = [0 im; -im 0]
2×2 Array{Complex{Int64},2}:
 0+0im  0+1im
 0-1im  0+0im

julia> eigmin(A)
-1.0

julia> A = [0 im; -1 0]
2×2 Array{Complex{Int64},2}:
  0+0im  0+1im
 -1+0im  0+0im

julia> eigmin(A)
ERROR: DomainError with Complex{Int64}[0+0im 0+1im; -1+0im 0+0im]:
`A` cannot have complex eigenvalues.
Stacktrace:
[...]
LinearAlgebra.eigvecs ── 関数
eigvecs(A::SymTridiagonal[, eigvals]) -> Matrix

A の固有ベクトルを列に持つ行列 M を返します。k 番目の固有ベクトルはスライス M[:, k] で取得できます。

省略可能な引数 eigvals に固有値のベクトルを指定すると、eigvecs はそれに対応する固有ベクトルを返します。

julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 1.0  2.0   
 2.0  2.0  3.0
     3.0  1.0

julia> eigvals(A)
3-element Array{Float64,1}:
 -2.1400549446402604
  1.0000000000000002
  5.140054944640259

julia> eigvecs(A)
3×3 Array{Float64,2}:
  0.418304  -0.83205      0.364299
 -0.656749  -7.39009e-16  0.754109
  0.627457   0.5547       0.546448

julia> eigvecs(A, [1.])
3×1 Array{Float64,2}:
  0.8320502943378438
  4.263514128092366e-17
 -0.5547001962252291
eigvecs(A; permute::Bool=true, scale::Bool=true, sortby) -> Matrix

A の固有ベクトルを列に持つ行列 M を返します。k 番目の固有ベクトルはスライス M[:, k] で取得できます。キーワード引数 permute, scale, sortbyeigen と同様です。

julia> eigvecs([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
eigvecs(A, B) -> Matrix

AB の一般化固有ベクトルを列に持つ行列 M を返します。k 番目の固有ベクトルはスライス M[:, k] で取得できます。

julia> A = [1 0; 0 -1]
2×2 Array{Int64,2}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Array{Int64,2}:
 0  1
 1  0

julia> eigvecs(A, B)
2×2 Array{Complex{Float64},2}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im
LinearAlgebra.eigen ── 関数
eigen(A; permute::Bool=true, scale::Bool=true, sortby) -> Eigen

A の固有値分解を計算します。返り値は分解を表す Eigen 型のオブジェクト F であり、固有値は F.values に、固有ベクトルは行列 F.vectors の列に含まれます (k 番目の固有ベクトルはスライス F.vectors[:, k] で取得できます)。

分解オブジェクトを反復すると F.valuesF.vectors が生成されます。

Eigen オブジェクトに対して利用可能な関数は inv, det, isposdef です。

対称でない一般の行列に対しては、固有値の計算を始める前に行列をどの程度調整するかを指定できます。オプション permute=true を指定すると上三角行列に近づけるために行列が置換され、scale=true を指定すると行列が対角要素でスケールされ行と列の絶対値が大まかに揃えられます。デフォルトでは両方とも true です。

デフォルトで固有値と固有ベクトルは (real(λ),imag(λ)) に従って辞書順にソートされますが、異なる比較関数 by(λ)sortby に渡すこともできます。また sortby=nothing とすれば固有値は任意の順序で並びます。DiagonalSymTridiagonal といった特殊な行列型は独自のソート規則を実装しており、キーワード引数 sortby を受け付けません。

julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
values:
3-element Array{Float64,1}:
  1.0
  3.0
 18.0
vectors:
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> F.values
3-element Array{Float64,1}:
  1.0
  3.0
 18.0

julia> F.vectors
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> vals, vecs = F; # 反復を使った分解の展開

julia> vals == F.values && vecs == F.vectors
true
eigen(A, B) -> GeneralizedEigen

AB の一般化固有値分解を計算します。返り値は分解を表す GeneralizedEigen 型のオブジェクト F であり、一般固有値は F.values に、一般固有ベクトルは行列 F.vectors の列に含まれます (k 番目の固有ベクトルはスライス F.vectors[:, k] で取得できます)。

分解オブジェクトを反復すると F.valuesF.vectors が生成されます。

eigen に渡されるキーワード引数は全て低水準の eigen! 関数に渡されます。

julia> A = [1 0; 0 -1]
2×2 Array{Int64,2}:
 1   0
 0  -1

julia> B = [0 1; 1 0]
2×2 Array{Int64,2}:
 0  1
 1  0

julia> F = eigen(A, B);

julia> F.values
2-element Array{Complex{Float64},1}:
 0.0 - 1.0im
 0.0 + 1.0im

julia> F.vectors
2×2 Array{Complex{Float64},2}:
  0.0+1.0im   0.0-1.0im
 -1.0+0.0im  -1.0-0.0im

julia> vals, vecs = F; # 反復を使った分解の展開

julia> vals == F.values && vecs == F.vectors
true
eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, irange::UnitRange) -> Eigen

A の固有値分解を計算します。返り値は分解を表す Eigen 型のオブジェクト F であり、固有値は F.values に、固有ベクトルは行列 F.vectors の列に含まれます (k 番目の固有ベクトルはスライス F.vectors[:, k] で取得できます)。

分解オブジェクトを反復すると F.valuesF.vectors が生成されます。

Eigen オブジェクトに対して利用可能な関数は inv, det, isposdef です。

UnitRange 型の値 irange は固有値をソートしたときの位置を表し、この値によって探索する固有値が指定されます。

情報

A の次元が n のとき irange1:n でないと、返り値の分解は途中で打ち切られた分解となります。

eigen(A::Union{SymTridiagonal, Hermitian, Symmetric}, vl::Real, vu::Real) -> Eigen

A の固有値分解を計算します。返り値は分解を表す Eigen 型のオブジェクト F であり、固有値は F.values に、固有ベクトルは行列 F.vectors の列に含まれます (k 番目の固有ベクトルはスライス F.vectors[:, k] で取得できます)。

分解オブジェクトを反復すると F.valuesF.vectors が生成されます。

Eigen オブジェクトに対して利用可能な関数は inv, det, isposdef です。

vl, vu は固有値を探索する範囲の下界および上界です。

情報

[vl, vu] に A の固有値が全て含まれない場合、返り値の分解は途中で打ち切られた分解となります

LinearAlgebra.eigen! ── 関数
eigen!(A, [B]; permute, scale, sortby)

基本的に eigen と同様ですが、入力 A, B のコピーを作らずに上書きすることで空間を節約します。

Hessenberg <: Factorization

正方行列のヘッセンベルグ分解 QHQ' またはそのシフト Q(H+μI)Q' を表します。hessenberg 関数がこの型の値を返します。

LinearAlgebra.hessenberg ── 関数
hessenberg(A) -> Hessenberg

A のヘッセンベルグ分解を計算し、Hessenberg オブジェクトを返します。F を分解オブジェクトとすれば、ユニタリ行列は F.Q (LinearAlgebra.HessenbergQ 型) であり、ヘッセンベルグ行列は F.H (UpperHessenberg 型) です。Matrix(F.H) あるいは Matrix(F.Q) とすれば通常の行列に変換できます。

A がエルミート (Hermitian) あるいは実対称 (Symmetric) のときヘッセンベルグ分解は実対称三重対角行列を生成し、F.HSymTridiagonal 型となります。

シフトされた分解 A+μI = Q (H+μI) Q'UniformScaling 型の I を使って F + μ*I として効率的に構築できることに注意してください。こうすると共有された格納場所と変更されたシフトを持つ新しい Hessenberg オブジェクトが作成されます。このため (F + μ*I) \ b という形のシフトされた求解が F を一度だけ構築することで高速に行えます。

分解オブジェクトを反復すると因子 F.Q, F.H, F.μ が生成されます。

julia> A = [4. 9. 7.; 4. 4. 1.; 4. 3. 2.]
3×3 Array{Float64,2}:
 4.0  9.0  7.0
 4.0  4.0  1.0
 4.0  3.0  2.0

julia> F = hessenberg(A)
Hessenberg{Float64,UpperHessenberg{Float64,Array{Float64,2}},Array{Float64,2},Array{Float64,1},Bool}
Q factor:
3×3 LinearAlgebra.HessenbergQ{Float64,Array{Float64,2},Array{Float64,1},false}:
 1.0   0.0        0.0
 0.0  -0.707107  -0.707107
 0.0  -0.707107   0.707107
H factor:
3×3 UpperHessenberg{Float64,Array{Float64,2}}:
  4.0      -11.3137       -1.41421
 -5.65685    5.0           2.0
           -8.88178e-16   1.0

julia> F.Q * F.H * F.Q'
3×3 Array{Float64,2}:
 4.0  9.0  7.0
 4.0  4.0  1.0
 4.0  3.0  2.0

julia> q, h = F; # 反復を使った分解の展開

julia> q == F.Q && h == F.H
true
hessenberg!(A) -> Hessenberg

hessenberg! は基本的に hessenberg と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。

LinearAlgebra.Schur ──
Schur <: Factorization

行列 A のシューア分解を表す型です。対応する行列分解関数 schur(A) はこの型の値を返します。

分解オブジェクト F::Schur の (疑) 三角のシューア因子は F.Shur または F.T で、直交/ユニタリのシューアベクトルは F.vectors または F.Z で取得でき、A = F.vectors*F.Schur*F.vectors' が成り立ちます。A の固有値は F.values で取得できます。

分解オブジェクトを反復すると F.T, F.Z, F.values が生成されます。

julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> F = schur(A)
Schur{Float64,Array{Float64,2}}
T factor:
2×2 Array{Float64,2}:
 3.0   9.0
 0.0  -2.0
Z factor:
2×2 Array{Float64,2}:
  0.961524  0.274721
 -0.274721  0.961524
eigenvalues:
2-element Array{Float64,1}:
  3.0
 -2.0

julia> F.vectors * F.Schur * F.vectors'
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> t, z, vals = F; # 反復を使った分解の展開

julia> t == F.T && z == F.Z && vals == F.values
true
GeneralizedSchur <: Factorization

二つの行列 AB の一般化シューア分解を表す型です。対応する行列分解関数 schur(A, B) はこの型の値を返します。

F::GeneralizedSchur を分解オブジェクトとすれば、(疑) 三角のシューア因子は F.SF.T で、左ユニタリ/直交シューアベクトルは F.left または F.Q で、右ユニタリ/直交シューアベクトルは F.right または F.Z で取得でき、A=F.left*F.S*F.right'B=F.left*F.T*F.right' が成り立ちます。AB の一般化固有値は F.α./F.β で取得できます。

分解オブジェクトを反復すると F.S, F.T, F.Q, F.Z, F.α, F.β が生成されます。

LinearAlgebra.schur ── 関数
schur(A::StridedMatrix) -> F::Schur

行列 A のシューア分解を返します。Shur 型のオブジェクト F の (疑) 三角のシューア因子は F.Schur あるいは F.T で、直交/ユニタリのシューアベクトルは F.vectors あるいは F.Z で取得でき、A = F.vectors*F.Schur*F.vectors' が成り立ちます。A の固有値は F.values で取得できます。

分解オブジェクトを反復すると F.T, F.Z, F.values が生成されます。

julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> F = schur(A)
Schur{Float64,Array{Float64,2}}
T factor:
2×2 Array{Float64,2}:
 3.0   9.0
 0.0  -2.0
Z factor:
2×2 Array{Float64,2}:
  0.961524  0.274721
 -0.274721  0.961524
eigenvalues:
2-element Array{Float64,1}:
  3.0
 -2.0

julia> F.vectors * F.Schur * F.vectors'
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> t, z, vals = F; # 反復を使った分解の展開

julia> t == F.T && z == F.Z && vals == F.values
true
schur(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur

行列 AB の一般化シューア分解 (QZ 分解) を計算します。FGeneralizedSchur オブジェクトとすれば、(疑) 三角のシューア因子は F.SF.T で、左ユニタリ/直交シューアベクトルは F.left または F.Q で、右ユニタリ/直交シューアベクトルは F.right または F.Z で取得でき、A=F.left*F.S*F.right'B=F.left*F.T*F.right' が成り立ちます。AB の一般化固有値は F.α./F.β で取得できます。

分解オブジェクトを反復すると F.S, F.T, F.Q, F.Z, F.α, F.β が生成されます。

LinearAlgebra.schur! ── 関数
schur!(A::StridedMatrix) -> F::Schur

基本的に schur と同じですが、引数 A を作業場所として使います。

julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0

julia> F = schur!(A)
Schur{Float64,Array{Float64,2}}
T factor:
2×2 Array{Float64,2}:
 3.0   9.0
 0.0  -2.0
Z factor:
2×2 Array{Float64,2}:
  0.961524  0.274721
 -0.274721  0.961524
eigenvalues:
2-element Array{Float64,1}:
  3.0
 -2.0

julia> A
2×2 Array{Float64,2}:
 3.0   9.0
 0.0  -2.0
schur!(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur

基本的に schur と同じですが、引数の行列 A, B を作業場所として使います。

LinearAlgebra.ordschur ── 関数
ordschur(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur

行列 A = Z*T*Z' のシューア分解 F を論理値配列 select に従って並べ替え、並べ替えた後の Schur オブジェクト F を返します。選択された固有値が F.Schur の対角要素の先頭部分となり、同じ個数の F.vectors の先頭の列が対応する左不変部分空間の直交/ユニタリな基底を構成します。A が実行列の場合には、固有値に含まれる複素共役の組は両方 select に含まれるか両方含まれないかのどちらかである必要があります。

ordschur(F::GeneralizedSchur,
         select::Union{Vector{Bool},BitVector}) -> F::GeneralizedSchur

行列の組 (A, B) = (Q*S*Z', Q*T*Z') の一般化シューア分解 F を受け取り、論理値配列 select に従って並べ替え、並べ替えた後の GeneralizedSchur オブジェクトを返します。選択された固有値が F.SF.T の対角要素の先頭部分となります。左および右の直交/ユニタリシューアベクトルは、並べ替えた後にも (A, B) = F.Q*(F.S, F.T)*F.Z' が成り立ち、さらに AB の一般化固有値が F.α./F.β で計算できるように並べ替えられます。

LinearAlgebra.ordschur! ── 関数
ordschur!(F::Schur, select::Union{Vector{Bool},BitVector}) -> F::Schur

基本的に ordschur と同じですが、分解 F に結果を上書きします。

ordschur!(F::GeneralizedSchur,
          select::Union{Vector{Bool},BitVector}) -> F::GeneralizedSchur

基本的に ordschur と同じですが、分解 F に結果を上書きします。

LinearAlgebra.SVD ──
SVD <: Factorization

行列 A の特異値分解 (SVD 分解) を表す型です。対応する行列分解関数 svd(A) はこの型の値を返します。

分解オブジェクト F::SVDU, S, V, Vt はそれぞれ F.U, F.S, F.V, F.Vt で取得でき、A = U * Diagonal(S) * Vt が成り立ちます。

分解オブジェクトを反復すると U, S, V が生成されます。

julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
4×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  2.0
 0.0  0.0  3.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  2.0  0.0  0.0  0.0

julia> F = svd(A)
SVD{Float64,Float64,Array{Float64,2}}
U factor:
4×4 Array{Float64,2}:
 0.0  1.0  0.0   0.0
 1.0  0.0  0.0   0.0
 0.0  0.0  0.0  -1.0
 0.0  0.0  1.0   0.0
singular values:
4-element Array{Float64,1}:
 3.0
 2.23606797749979
 2.0
 0.0
Vt factor:
4×5 Array{Float64,2}:
 -0.0       0.0  1.0  -0.0  0.0
  0.447214  0.0  0.0   0.0  0.894427
 -0.0       1.0  0.0  -0.0  0.0
  0.0       0.0  0.0   1.0  0.0

julia> F.U * Diagonal(F.S) * F.Vt
4×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  2.0
 0.0  0.0  3.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  2.0  0.0  0.0  0.0

julia> u, s, v = F; # 反復を使った分解の展開

julia> u == F.U && s == F.S && v == F.V
true
GeneralizedSVD <: Factorization

二つの行列 A, B の一般化特異値分解 (一般化 SVD 分解) を表す型です。A = F.U*F.D1*F.R0*F.Q'B = F.V*F.D2*F.R0*F.Q' が成り立ちます。対応する行列分解関数 svd(A, B) はこの型の値を返します。

M×N 行列 AP×N 行列 B に対して:

  • UM×M の直交行列です。
  • VP×P の直交行列です。
  • QN×N の直交行列です。
  • D1M×(K+L) の直交行列で、対角要素の最初の K 個は 1 です。
  • D2P×(K+L) 行列で、右上の L×L のブロックが対角行列です。
  • R0(K+L)×N 行列で、右側の (K+L)×(K+L) ブロックが非特異のブロック上三角行列です。

K+L は行列 [A; B] の実効数値ランク (effective numerical rank) です。

分解オブジェクトを反復すると U, V, Q, D1, D2, R0 が生成されます。

LAPACK のユーザーガイドで説明されるように、F.D1F.D2 には関連性があります。この関数は xGGSVD3 ルーチンを内部で呼び出します (LAPACK 3.6.0 以降)。

julia> A = [1. 0.; 0. -1.]
2×2 Array{Float64,2}:
 1.0   0.0
 0.0  -1.0

julia> B = [0. 1.; 1. 0.]
2×2 Array{Float64,2}:
 0.0  1.0
 1.0  0.0

julia> F = svd(A, B)
GeneralizedSVD{Float64,Array{Float64,2}}
U factor:
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
V factor:
2×2 Array{Float64,2}:
 -0.0  -1.0
  1.0   0.0
Q factor:
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
D1 factor:
2×2 SparseArrays.SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 1]  =  0.707107
  [2, 2]  =  0.707107
D2 factor:
2×2 SparseArrays.SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 1]  =  0.707107
  [2, 2]  =  0.707107
R0 factor:
2×2 Array{Float64,2}:
 1.41421   0.0
 0.0      -1.41421

julia> F.U*F.D1*F.R0*F.Q'
2×2 Array{Float64,2}:
 1.0   0.0
 0.0  -1.0

julia> F.V*F.D2*F.R0*F.Q'
2×2 Array{Float64,2}:
 0.0  1.0
 1.0  0.0
LinearAlgebra.svd ── 関数
svd(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD

A の特異値分解を計算し、SVD オブジェクトを返します。

分解オブジェクト F から U, S, V, Vt を得るには、それぞれ F.U, F.S, F.V, F.Vt とします。このとき A = U*Diagonal(S)*Vt が成り立ちます。使われるアルゴリズムは Vt を生成するので、V を取得するよりも Vt を取得した方が効率的です。S に含まれる特異値は降順にソートされます。

分解オブジェクトを反復すると U, S, V が生成されます。

full = false (デフォルト) だと “薄い” SVD 分解が返ります。完全な分解では UM×M 行列で VN×N 行列となりますが、薄い分解では UM×K 行列で VN×K 行列となります。ここで K = min(M,N) であり、特異値の個数を表します。

alg = DivideAndConquer() だと分割統治アルゴリズムを使って SVD が計算されます。もう一つの選択肢は alg = QRIteration() です (通常こちらの方が低速ですが、高精度になります)。

Julia 1.3

キーワード引数 alg は Julia 1.3 以降でサポートされます。

julia> A = rand(4,3);

julia> F = svd(A); # 分解を表すオブジェクトの保存

julia> A  F.U * Diagonal(F.S) * F.Vt
true

julia> U, S, V = F; # 反復を使った分解の展開

julia> A  U * Diagonal(S) * V'
true

julia> Uonly, = svd(A); # U だけを格納

julia> Uonly == U
true
svd(A, B) -> GeneralizedSVD

AB の一般化 SVD を計算し、GeneralizedSVD オブジェクト F を返します。返り値に対して [A;B] = [F.U * F.D1; F.V * F.D2] * F.R0 * F.Q' が成り立ちます。

  • UM×M の直交行列です。
  • VP×P の直交行列です。
  • QN×N の直交行列です。
  • D1M×(K+L) の直交行列で、対角要素の最初の K 個は 1 です。
  • D2P×(K+L) 行列で、右上の L×L のブロックが対角行列です。
  • R0(K+L)×N 行列で、右側の (K+L)×(K+L) ブロックが非特異のブロック上三角行列です。

K+L は行列 [A; B] の実効数値ランク (effective numerical rank) です。

分解オブジェクトを反復すると U, V, Q, D1, D2, R0 が生成されます。

一般化 SVD は現実の問題において A に属する量と B に属する量を比較するような場合に使われます。例えば人間の遺伝子とイーストの遺伝子、信号と雑音、クラスター間とクラスター内などの比較で使われます (詳しい議論は Edelman と Wang による The GSVD: Where are the ellipses? を参照してください)。

一般化 SVD は [A; B][UC; VS]H に分解します。ここで [UC; VS][A; B] の列空間に対する自然な直交基底、H = RQ'[A;B] の行空間に対する自然な非直交基底であり、H の最も上にある行が最も強く行列 A に結び付き、最も下にある行が最も強く行列 B に結び付きます。マルチコサイン/サイン行列 C, SA の量と B の量のマルチメジャー (multi-measure) を表し、UV はそのマルチメジャーが測られた方向を表します。

julia> A = randn(3,2); B=randn(4,2);

julia> F = svd(A, B);

julia> U,V,Q,C,S,R = F;

julia> H = R*Q';

julia> [A; B]  [U*C; V*S]*H
true

julia> [A; B]  [F.U*F.D1; F.V*F.D2]*F.R0*F.Q'
true

julia> Uonly, = svd(A,B);

julia> U == Uonly
true
LinearAlgebra.svd! ── 関数
svd!(A; full::Bool = false, alg::Algorithm = default_svd_alg(A)) -> SVD

基本的に svd と同じですが、入力 A のコピーを作らずに上書きすることで空間を節約します。

svd!(A, B) -> GeneralizedSVD

基本的に svd と同様ですが、引数 AB のコピーを作らずにインプレースに改変します。詳細は svd のドキュメントを見てください。

LinearAlgebra.svdvals ── 関数
svdvals(A)

A の固有値を降順に返します。

julia> A = [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
4×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  2.0
 0.0  0.0  3.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  2.0  0.0  0.0  0.0

julia> svdvals(A)
4-element Array{Float64,1}:
 3.0
 2.23606797749979
 2.0
 0.0
svdvals(A, B)

AB の一般化固有値分解から一般化固有値を求めます。svd も参照してください。

julia> A = [1. 0.; 0. -1.]
2×2 Array{Float64,2}:
 1.0   0.0
 0.0  -1.0

julia> B = [0. 1.; 1. 0.]
2×2 Array{Float64,2}:
 0.0  1.0
 1.0  0.0

julia> svdvals(A, B)
2-element Array{Float64,1}:
 1.0
 1.0
LinearAlgebra.svdvals! ── 関数
svdvals!(A)

A を上書きすることで空間を節約しながら A の固有値を求めます。svdvalssvd も参照してください。

svdvals!(A, B)

AB の一般化固有値分解から一般化固有値を求めます。そのとき AB を上書きすることで空間を節約します。svdsvdvals も参照してください。

LinearAlgebra.Givens(i1,i2,c,s) -> G

ギブンス回転を表す線形作用素です。cs はそれぞれ回転角のコサインとサインを表します。Givens 型は左乗算 G*A と共役転置の右乗算 A*G' をサポートします。この型は size を持たないので、乗算が意味を持つサイズを持つ任意の行列との乗算が可能です。具体的には i2<=size(A,2) を満たす全ての A に対して G*A が、i2<=size(A,1) を満たす全ての A に対して A*G' が行えます。

givens も参照してください。

LinearAlgebra.givens ── 関数
givens(f::T, g::T, i1::Integer, i2::Integer) where {T} -> (G::Givens, r::T)

ギブンス回転 G とスカラー r を計算します。返り値は次の条件を満たします:

x[i1] = f
x[i2] = g

が成り立つ任意の x に対して、乗算

y = G*x

で定義される y が次の性質を持つ:

y[i1] = r
y[i2] = 0

LinearAlgebra.Givens も参照してください。

givens(A::AbstractArray, i1::Integer, i2::Integer, j::Integer) -> (G::Givens, r)

ギブンス回転 G とスカラー r を計算します。返り値は次の条件を満たします: 積

B = G*A

が、次の性質を持つ:

B[i1,j] = r
B[i2,j] = 0

LinearAlgebra.Givens も参照してください。

givens(x::AbstractVector, i1::Integer, i2::Integer) -> (G::Givens, r)

ギブンス回転 G とスカラー r を計算します。返り値は次の条件を満たします: 積

B = G*x

が、次の性質を持つ:

B[i1] = r
B[i2] = 0

LinearAlgebra.Givens も参照してください。

LinearAlgebra.triu ── 関数
triu(M)

行列 M の上三角部分からなる上三角行列を返します。

julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> triu(a)
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0
 0.0  0.0  0.0  1.0
triu(M, k::Integer)

M の第 k 優対角要素およびそれより上の部分からなる上三角行列を返します。

julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> triu(a,3)
4×4 Array{Float64,2}:
 0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> triu(a,-3)
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
LinearAlgebra.triu! ── 関数
triu!(M)

M を上三角行列に上書きして返します。triu も参照してください。

triu!(M, k::Integer)

M を第 k 優対角要素およびそれより上の部分からなる上三角行列に上書きして返します。

julia> M = [1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5]
5×5 Array{Int64,2}:
 1  2  3  4  5
 1  2  3  4  5
 1  2  3  4  5
 1  2  3  4  5
 1  2  3  4  5

julia> triu!(M, 1)
5×5 Array{Int64,2}:
 0  2  3  4  5
 0  0  3  4  5
 0  0  0  4  5
 0  0  0  0  5
 0  0  0  0  0
LinearAlgebra.tril ── 関数
tril(M)

行列 M の下三角部分からなる行列を返します。

julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> tril(a)
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 1.0  1.0  0.0  0.0
 1.0  1.0  1.0  0.0
 1.0  1.0  1.0  1.0
tril(M, k::Integer)

M の第 k 優対角要素およびそれより下の部分からなる下三角行列を返します。

julia> a = fill(1.0, (4,4))
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> tril(a,3)
4×4 Array{Float64,2}:
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0

julia> tril(a,-3)
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0
LinearAlgebra.tril! ── 関数
tril!(M)

M を下三角行列に上書きして返します。tril も参照してください。

tril!(M, k::Integer)

M を第 k 優対角要素およびそれより下の部分からなる下三角行列に上書きして返します。

julia> M = [1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5; 1 2 3 4 5]
5×5 Array{Int64,2}:
 1  2  3  4  5
 1  2  3  4  5
 1  2  3  4  5
 1  2  3  4  5
 1  2  3  4  5

julia> tril!(M, 2)
5×5 Array{Int64,2}:
 1  2  3  0  0
 1  2  3  4  0
 1  2  3  4  5
 1  2  3  4  5
 1  2  3  4  5
LinearAlgebra.diagind ── 関数
diagind(M, k::Integer=0)

行列 M の第 k 優対角要素の添え字を与える AbstractRange を返します。

julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> diagind(A,-1)
2:4:6
LinearAlgebra.diag ── 関数
diag(M, k::Integer=0)

行列 M の第 k 優対角要素をベクトルとして返します。

diagm も参照してください。

julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> diag(A,1)
2-element Array{Int64,1}:
 2
 6
LinearAlgebra.diagm ── 関数
diagm(kv::Pair{<:Integer,<:AbstractVector}...)
diagm(m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)

整数とベクトルの Pair から行列を構築します。ベクトル kv.second が第 kv.first 優対角要素となります。デフォルトでは返り値の行列は正方であり、そのサイズは kv から計算されます。正方でないサイズ m×n は引数の最初に m,n を渡すことで指定できます (空いた部分はゼロで埋められます)。

diagm は完全な行列を構築します。空間効率に優れ高速な算術が可能なバージョンが必要なときは Diagonal, BidiagonalTridiagonal, SymTridiagonal を使ってください。

julia> diagm(1 => [1,2,3])
4×4 Array{Int64,2}:
 0  1  0  0
 0  0  2  0
 0  0  0  3
 0  0  0  0

julia> diagm(1 => [1,2,3], -1 => [4,5])
4×4 Array{Int64,2}:
 0  1  0  0
 4  0  2  0
 0  5  0  3
 0  0  0  0
diagm(v::AbstractVector)
diagm(m::Integer, n::Integer, v::AbstractVector)

ベクトルの要素を対角要素に持つ行列を構築します。デフォルト (size = nothing) では返り値の行列は正方であり、そのサイズは length(v) から定まります。正方でないサイズ m×n を持つ行列は引数の最初に m,n を渡すことで指定できます。

julia> diagm([1,2,3])
3×3 Array{Int64,2}:
 1  0  0
 0  2  0
 0  0  3
LinearAlgebra.rank ── 関数
rank(A::AbstractMatrix; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
rank(A::AbstractMatrix, rtol::Real)

行列 A のランクを計算します。ランクの計算は A の特異値で max(atol, rtol*σ₁) より大きいものの個数を数えることで行われます。ここで σ₁A の最大特異値で、atolrtol はそれぞれ絶対的および相対的な許容値です。デフォルトの相対的許容値は n*ϵ であり、ここで nA の小さい方の次元のサイズ、ϵA の要素型の eps を表します。

Julia 1.1

キーワード引数 atol, rtol は Julia 1.1 以降でサポートされます。Julia 1.0 では rtol が位置引数として利用可能でしたが、これは Julia 2.0 で非推奨となります。

julia> rank(Matrix(I, 3, 3))
3

julia> rank(diagm(0 => [1, 0, 2]))
2

julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.1)
2

julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.00001)
3

julia> rank(diagm(0 => [1, 0.001, 2]), atol=1.5)
1
LinearAlgebra.norm ── 関数
norm(A, p::Real=2)

数値 (あるいは norm が定義される) 型の値を生成する反復可能コンテナ A を受け取り、A をコンテナと同じ長さのベクトルとみなしたときの p ノルムを計算します。デフォルトでは p=2 で、A は任意次元の配列でも構いません。

p ノルムは次のように定義されます:

\[ \|A\|_p = \left( \sum_{i=1}^n | a_i | ^p \right)^{1/p} \]

ここで \(a_i\)\(A\) の要素、\(|a_i|\)\(a_i\)norm です。この p ノルムは A の要素に対する (p=2 の) norm を使って計算されるので、p != 2 のときベクトルのベクトルの p ノルムはそれをブロックベクトルとして扱ったときの p ノルムと一般に互換性がありません。

p は任意の数値にできます (ただし全ての値が数学的に正当なベクトルノルムを生成するとは限りません)。特に norm(A, Inf)abs.(A) の最大値を返し、norm(A, -Inf)abs.(A) の最小値を返します。A が行列で p=2 なら、この関数の返り値はフロベニウスノルムと一致します。

第二引数 pnorm に対するインターフェースの一部になっていなくても構いません。つまり、独自型が第二引数を持たない norm(A) だけを実装することは許されます。

行列の作用素ノルムを計算するには opnorm を使ってください。

julia> v = [3, -2, 6]
3-element Array{Int64,1}:
  3
 -2
  6

julia> norm(v)
7.0

julia> norm(v, 1)
11.0

julia> norm(v, Inf)
6.0

julia> norm([1 2 3; 4 5 6; 7 8 9])
16.881943016134134

julia> norm([1 2 3 4 5 6 7 8 9])
16.881943016134134

julia> norm(1:9)
16.881943016134134

julia> norm(hcat(v,v), 1) == norm(vcat(v,v), 1) != norm([v,v], 1)
true

julia> norm(hcat(v,v), 2) == norm(vcat(v,v), 2) == norm([v,v], 2)
true

julia> norm(hcat(v,v), Inf) == norm(vcat(v,v), Inf) != norm([v,v], Inf)
true
norm(x::Number, p::Real=2)

数値に対する norm\(\left( |x|^p \right)^{1/p}\) を返します。

julia> norm(2, 1)
2.0

julia> norm(-2, 1)
2.0

julia> norm(2, 2)
2.0

julia> norm(-2, 2)
2.0

julia> norm(2, Inf)
2.0

julia> norm(-2, Inf)
2.0
LinearAlgebra.opnorm ── 関数
opnorm(A::AbstractMatrix, p::Real=2)

p ノルムによって誘導される作用素ノルム (行列ノルム) を返します。p として正当な値は 1, 2, Inf のみです (注意: 疎行列に対する p=2 のノルムは現在実装されていません)。フロベニウスノルムの計算には norm を使ってください。

p=1 のとき、作用素ノルムは A の列ごとに計算した絶対値の和の最大値を返します:

\[ \|A\|_1 = \max_{1 ≤ j ≤ n} \sum_{i=1}^m | a_{ij} | \]

ここで \(a_{ij}\)\(A\) の要素で、\(m\)\(n\) はその次元です。

p=2 の作用素ノルムはスペクトルノルムであり、A の最大特異値と一致します。

p=Inf の作用素ノルムは A の行ごとに計算した絶対値の和の最大値です:

\[ \|A\|_\infty = \max_{1 ≤ i ≤ m} \sum _{j=1}^n | a_{ij} | \]

julia> A = [1 -2 -3; 2 3 -1]
2×3 Array{Int64,2}:
 1  -2  -3
 2   3  -1

julia> opnorm(A, Inf)
6.0

julia> opnorm(A, 1)
5.0
opnorm(x::Number, p::Real=2)

数値に対する opnorm\(\left( |x|^p \right)^{1/p}\) を返します。これは norm と等価です。

opnorm(A::Adjoint{<:Any,<:AbstracVector}, q::Real=2)
opnorm(A::Transpose{<:Any,<:AbstracVector}, q::Real=2)

随伴または転置でラップされたベクトル A に対する opnormAq ノルムを返します。これは p = q/(q-1) に対する p ノルムと等価です。q = 2 のとき p = q = 2 となって二つのノルムは一致します。ベクトルとしての A に対する p ノルムの計算では norm が使われます。

ベクトル空間とその双対空間におけるノルムの違いは、双対性とドット積の関係を保存するために存在します。返り値は 1×n 行列の作用素 p ノルムと一致します。

julia> v = [1; im];

julia> vc = v';

julia> opnorm(vc, 1)
1.0

julia> norm(vc, 1)
2.0

julia> norm(v, 1)
2.0

julia> opnorm(vc, 2)
1.4142135623730951

julia> norm(vc, 2)
1.4142135623730951

julia> norm(v, 2)
1.4142135623730951

julia> opnorm(vc, Inf)
2.0

julia> norm(vc, Inf)
1.0

julia> norm(v, Inf)
1.0
LinearAlgebra.normalize! ── 関数
normalize!(a::AbstractArray, p::Real=2)

配列 a をインプレースに正規化して、p ノルムが 1 に等しいように、つまり norm(a, p) == 1 となるようにします。normalizenorm も参照してください。

LinearAlgebra.normalize ── 関数
normalize(a::AbstractArray, p::Real=2)

配列 a を正規化して、p ノルムが 1 に等しいように、つまり norm(a, p) == 1 となるようにしたものを返します。normalize!norm も参照してください。

julia> a = [1,2,4];

julia> b = normalize(a)
3-element Array{Float64,1}:
 0.2182178902359924
 0.4364357804719848
 0.8728715609439696

julia> norm(b)
1.0

julia> c = normalize(a, 1)
3-element Array{Float64,1}:
 0.14285714285714285
 0.2857142857142857
 0.5714285714285714

julia> norm(c, 1)
1.0

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

julia> norm(a)
6.48074069840786

julia> normalize(a)
2×3 Array{Float64,2}:
 0.154303  0.308607  0.617213
 0.154303  0.308607  0.617213

LinearAlgebra.cond ── 関数
cond(M, p::Real=2)

作用素 p ノルムを使って計算した行列 M の条件数を返します。p に指定できるのは 1, 2, Inf であり、2 がデフォルトです。

LinearAlgebra.condskeel ── 関数
condskeel(M, [x, p::Real=Inf])
\[ \kappa_S(M, p) = \left\Vert \left\vert M \right\vert \left\vert M^{-1} \right\vert \right\Vert_p \\[10pt] \kappa_S(M, x, p) = \frac{\left\Vert \left\vert M \right\vert \left\vert M^{-1} \right\vert \left\vert x \right\vert \right\Vert_p}{\left \Vert x \right \Vert_p} \]

作用素 p ノルムを使って計算した行列 M の Skeel 条件数 \(\kappa_S\) を計算します。ベクトル x を与えれば、それに関する Skeel 条件数が計算されます。\(\left\vert M \right\vert\)M の各要素をその絶対値に置き換えた行列を表します: つまり \(\left\vert M \right\vert_{ij} = \left\vert M_{ij} \right\vert\) です。p に指定できるのは 1, 2, Inf であり、2 がデフォルトです。

この値は文献によっては Bauer 条件数・相対 (relative) 条件数・要素単位相対 (componentwise relative) 条件数とも呼ばれます。

LinearAlgebra.tr ── 関数
tr(M)

行列のトレースです。M の対角要素の和が返ります。

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

julia> tr(A)
5
LinearAlgebra.det ── 関数
det(M)

行列式を計算します。

julia> M = [1 0; 2 2]
2×2 Array{Int64,2}:
 1  0
 2  2

julia> det(M)
2.0
LinearAlgebra.logdet ── 関数
logdet(M)

行列式の対数です。log(det(M)) と等価ですが、精度や速度が向上する可能性があります。

julia> M = [1 0; 2 2]
2×2 Array{Int64,2}:
 1  0
 2  2

julia> logdet(M)
0.6931471805599453

julia> logdet(Matrix(I, 3, 3))
0.0
LinearAlgebra.logabsdet ── 関数
logabsdet(M)

行列式の絶対値の対数です。(log(abs(det(M))), sign(det(M))) と等価ですが、精度や速度が向上する可能性があります。

julia> A = [-1. 0.; 0. 1.]
2×2 Array{Float64,2}:
 -1.0  0.0
  0.0  1.0

julia> det(A)
-1.0

julia> logabsdet(A)
(0.0, -1.0)

julia> B = [2. 0.; 0. 1.]
2×2 Array{Float64,2}:
 2.0  0.0
 0.0  1.0

julia> det(B)
2.0

julia> logabsdet(B)
(0.6931471805599453, 1.0)
Base.inv ── メソッド
inv(M)

逆行列を返します。恒等行列 I に対して M * N = I が成り立つような行列 N を計算します。計算は左除算 N = M \ I を解くことで行われます。

julia> M = [2 5; 1 3]
2×2 Array{Int64,2}:
 2  5
 1  3

julia> N = inv(M)
2×2 Array{Float64,2}:
  3.0  -5.0
 -1.0   2.0

julia> M*N == N*M == Matrix(I, 2, 2)
true
LinearAlgebra.pinv ── 関数
pinv(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
pinv(M, rtol::Real) = pinv(M; rtol=rtol) # Julia 2.0 で非推奨

ムーア-ペンローズ疑似逆行列を計算します。

浮動小数点数を要素に持つ行列 M に対しては、max(atol, rtol*σ₁) より大きい特異値だけを逆転させる疑似逆行列の計算が有用です。ここで σ₁M の最大特異値を表します。

絶対的な許容値 (atol) および相対的な許容値 (rtol) に対する最適な値は M の値と疑似逆行列の用途の両方に応じて異なります。デフォルトの相対的な許容値は n*ϵ です。ここで nM の一番小さな次元のサイズ、ϵM の要素型の eps を表します。

悪条件の行列に対して最小二乗の意味で逆数を取るときは

rtol = sqrt(eps(real(float(one(eltype(M))))))

が推奨されます。

詳しくは [issue8859], [B96], [S84], [KY88] を参照してください。

julia> M = [1.5 1.3; 1.2 1.9]
2×2 Array{Float64,2}:
 1.5  1.3
 1.2  1.9

julia> N = pinv(M)
2×2 Array{Float64,2}:
  1.47287   -1.00775
 -0.930233   1.16279

julia> M * N
2×2 Array{Float64,2}:
 1.0          -2.22045e-16
 4.44089e-16   1.0
LinearAlgebra.nullspace ── 関数
nullspace(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
nullspace(M, rtol::Real) = nullspace(M; rtol=rtol) # Julia 2.0 で非推奨

M の零空間の基底を、絶対値が max(atol, rtol*σ₁) より大きい特異値に対応する M の特異ベクトルを考えることで計算します。ここで σ₁M の最大特異値です。

デフォルトの相対的な許容値は n*ϵ です。ここで nM の一番小さな次元のサイズ、ϵM の要素型の eps を表します。

julia> M = [1 0 0; 0 1 0; 0 0 0]
3×3 Array{Int64,2}:
 1  0  0
 0  1  0
 0  0  0

julia> nullspace(M)
3×1 Array{Float64,2}:
 0.0
 0.0
 1.0

julia> nullspace(M, rtol=3)
3×3 Array{Float64,2}:
 0.0  1.0  0.0
 1.0  0.0  0.0
 0.0  0.0  1.0

julia> nullspace(M, atol=0.95)
3×1 Array{Float64,2}:
 0.0
 0.0
 1.0
Base.kron ── 関数
kron(A, B)

二つのベクトルまたは二つの行列のクロネッカーテンソル積です。

実ベクトル vw に対して、クロネッカー積と直積には kron(v,w) == vec(w * transpose(v)) および w * transpose(v) == reshape(kron(v,w), (length(w), length(v))) という関係があります。vw の順序が等式の左右で入れ替わっているのに注意してください (データが列優先の順序で格納されるためです)。複素ベクトルに対する直積 w * v'v の共役を取る点でクロネッカー積 reshape(kron(v,w), (length(w), length(v))) と異なります。

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

julia> B = [im 1; 1 -im]
2×2 Array{Complex{Int64},2}:
 0+1im  1+0im
 1+0im  0-1im

julia> kron(A, B)
4×4 Array{Complex{Int64},2}:
 0+1im  1+0im  0+2im  2+0im
 1+0im  0-1im  2+0im  0-2im
 0+3im  3+0im  0+4im  4+0im
 3+0im  0-3im  4+0im  0-4im

julia> v = [1, 2]; w = [3, 4, 5];

julia> w*transpose(v)
3×2 Array{Int64,2}:
 3   6
 4   8
 5  10

julia> reshape(kron(v,w), (length(w), length(v)))
3×2 Array{Int64,2}:
 3   6
 4   8
 5  10
Base.exp ── メソッド
exp(A::AbstractMatrix)

A の行列指数を計算します。次のように定義されます:

\[ e^A = \sum_{n=0}^{\infty} \frac{A^n}{n!} \]

対称/エルミート行列 A に対しては固有値分解 (eigen) が使われます。それ以外の場合はスケーリングと二乗を用いたアルゴリズム [H05] が使われます。

julia> A = Matrix(1.0I, 2, 2)
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

julia> exp(A)
2×2 Array{Float64,2}:
 2.71828  0.0
 0.0      2.71828
Base.:^ ── メソッド
^(A::AbstractMatrix, p::Number)

行列のべき乗です。\(\exp(p\log(A))\) と等価です。

julia> [1 2; 0 3]^3
2×2 Array{Int64,2}:
 1  26
 0  27
Base.:^ ── メソッド
^(b::Number, A::AbstractMatrix)

行列指数関数です。\(\exp(\log(b)A)\) と等価です。

Julia 1.1

Irrational (例えば ) の行列指数は Julia 1.1 以降でサポートされます。

julia> 2^[1 2; 0 3]
2×2 Array{Float64,2}:
 2.0  6.0
 0.0  8.0

julia> ^[1 2; 0 3]
2×2 Array{Float64,2}:
 2.71828  17.3673
 0.0      20.0855
Base.log ── メソッド
log(A{T}::StridedMatrix{T})

A が負の実固有値を持たないなら、A の行列対数の主値を計算します。つまり \(e^X = A\) かつ \(X\) の全ての固有値 \(\lambda\) に対して \(-\pi < \text{Im}(\lambda) < \pi\) が成り立つような唯一の \(X\) が返ります。A が正でない実固有値を持つときは、可能なら主値でない行列対数が返ります。

対称/エルミート行列では固有値分解 (eigen) が使われ、三角行列では逆スケーリングと二乗を用いたアルゴリズム [AHR13] が使われます。一般の行列に対しては複素シューア分解 (schur) を計算してから三角因子に三角行列のアルゴリズムが適用されます。

julia> A = Matrix(2.7182818*I, 2, 2)
2×2 Array{Float64,2}:
 2.71828  0.0
 0.0      2.71828

julia> log(A)
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
Base.sqrt ── メソッド
sqrt(A::AbstractMatrix)

A が負の実固有値を持たないなら、A の行列平方根の主値を計算します。つまり \(X^2 = A\) を満たし固有値の実部が正であるような唯一の \(X\) が返ります。そうでなければ主値でない行列平方根が返ります。

A が実対称またはエルミートなら、固有値分解 (eigen) を使って平方根が計算されます。そういった行列では丸め誤差により生じるごく小さな負の固有値 λ はゼロとして扱われます。正確に言うと、全ての固有値が -rtol*(max |λ|) 以上の行列は半正定値で負の固有値は 0 であるとみなされます (このとき平方根はエルミートです)。 rtolA が実対称/エルミートのときにだけ sqrt に指定できるキーワード引数であり、デフォルトではマシン精度に size(A,1) を乗じた値が使われます。

それ以外の場合には、平方根は Björck-Hammarling 法 [BH83] を使って計算されます。これは複素シューア形式 (schur) を計算し、三角因子の平方根を取る方法です。

julia> A = [4 0; 0 4]
2×2 Array{Int64,2}:
 4  0
 0  4

julia> sqrt(A)
2×2 Array{Float64,2}:
 2.0  0.0
 0.0  2.0
Base.cos ── メソッド
cos(A::AbstractMatrix)

正方行列 A の行列コサインを計算します。

行列コサインは A が対称またはエルミートなら固有値分解 (eigen) を使って計算され、それ以外の場合には exp を呼び出すことで計算されます。

julia> cos(fill(1.0, (2,2)))
2×2 Array{Float64,2}:
  0.291927  -0.708073
 -0.708073   0.291927
Base.sin ── メソッド
sin(A::AbstractMatrix)

正方行列 A の行列サインを計算します。

行列サインは A が対称またはエルミートなら固有値分解 (eigen) を使って計算され、それ以外の場合には exp を呼び出すことで計算されます。

julia> sin(fill(1.0, (2,2)))
2×2 Array{Float64,2}:
 0.454649  0.454649
 0.454649  0.454649
Base.Math.sincos ── メソッド
sincos(A::AbstractMatrix)

正方行列 A の行列サインと行列コサインを計算します。

julia> S, C = sincos(fill(1.0, (2,2)));

julia> S
2×2 Array{Float64,2}:
 0.454649  0.454649
 0.454649  0.454649

julia> C
2×2 Array{Float64,2}:
  0.291927  -0.708073
 -0.708073   0.291927
Base.tan ── メソッド
tan(A::AbstractMatrix)

正方行列 A の行列タンジェントを計算します。

行列タンジェントは A が対称またはエルミートなら固有値分解 (eigen) を使って計算され、それ以外の場合には exp を呼び出すことで計算されます。

julia> tan(fill(1.0, (2,2)))
2×2 Array{Float64,2}:
 -1.09252  -1.09252
 -1.09252  -1.09252
Base.Math.sec ── メソッド
sec(A::AbstractMatrix)

正方行列 A の行列セカントを計算します。

Base.Math.csc ── メソッド
csc(A::AbstractMatrix)

正方行列 A の行列コセカントを計算します。

Base.Math.cot ── メソッド
cot(A::AbstractMatrix)

正方行列 A の行列コタンジェントを計算します。

Base.cosh ── メソッド
cosh(A::AbstractMatrix)

正方行列 A の行列ハイパボリックコサインを計算します。

Base.sinh ── メソッド
sinh(A::AbstractMatrix)

正方行列 A の行列ハイパボリックサインを計算します。

Base.tanh ── メソッド
tanh(A::AbstractMatrix)

正方行列 A の行列ハイパボリックタンジェントを計算します。

Base.Math.sech ── メソッド
sech(A::AbstractMatrix)

正方行列 A の行列ハイパボリックセカントを計算します。

Base.Math.csch ── メソッド
csch(A::AbstractMatrix)

正方行列 A の行列ハイパボリックコセカントを計算します。

Base.Math.coth ── メソッド
coth(A::AbstractMatrix)

正方行列 A の行列ハイパボリックコタンジェントを計算します。

Base.acos ── メソッド
acos(A::AbstractMatrix)

行列コサインの逆関数を正方行列 A に適用します。

逆コサインは A が対称またはエルミートなら固有値分解 (eigen) を使って計算され、それ以外の場合には logsqrt を呼び出すことで計算されます。この関数の計算で使われる理論と対数公式については [AH16] を参照してください。

julia> acos(cos([0.5 0.1; -0.2 0.3]))
2×2 Array{Complex{Float64},2}:
  0.5-8.32667e-17im  0.1+0.0im
 -0.2+2.63678e-16im  0.3-3.46945e-16im
Base.asin ── メソッド
asin(A::AbstractMatrix)

行列サインの逆関数を正方行列 A に適用します。

逆サインは A が対称またはエルミートなら固有値分解 (eigen) を使って計算され、それ以外の場合には logsqrt を呼び出すことで計算されます。この関数の計算で使われる理論と対数公式については [AH16] を参照してください。

julia> asin(sin([0.5 0.1; -0.2 0.3]))
2×2 Array{Complex{Float64},2}:
  0.5-4.16334e-17im  0.1-5.55112e-17im
 -0.2+9.71445e-17im  0.3-1.249e-16im
Base.atan ── メソッド
atan(A::AbstractMatrix)

行列タンジェントの逆関数を正方行列 A に適用します。

逆タンジェントは A が対称またはエルミートなら固有値分解 (eigen) を使って計算され、それ以外の場合には logsqrt を呼び出すことで計算されます。この関数の計算で使われる理論と対数公式については [AH16] を参照してください。

julia> atan(tan([0.5 0.1; -0.2 0.3]))
2×2 Array{Complex{Float64},2}:
  0.5+1.38778e-17im  0.1-2.77556e-17im
 -0.2+6.93889e-17im  0.3-4.16334e-17im
Base.Math.asec ── メソッド
asec(A::AbstractMatrix)

行列セカントの逆関数を正方行列 A に適用します。

Base.Math.acsc ── メソッド
acsc(A::AbstractMatrix)

行列コセカントの逆関数を正方行列 A に適用します。

Base.Math.acot ── メソッド
acot(A::AbstractMatrix)

行列コタンジェントの逆関数を正方行列 A に適用します。

Base.acosh ── メソッド
acosh(A::AbstractMatrix)

行列ハイパボリックコサインの逆関数を正方行列 A に適用します。この関数の計算で使われる理論と対数公式については [AH16] を参照してください。

Base.asinh ── メソッド
asinh(A::AbstractMatrix)

行列ハイパボリックサインの逆関数を正方行列 A に適用します。この関数の計算で使われる理論と対数公式については [AH16] を参照してください。

Base.atanh ── メソッド
atanh(A::AbstractMatrix)

行列ハイパボリックタンジェントの逆関数を正方行列 A に適用します。この関数の計算で使われる理論と対数公式については [AH16] を参照してください。

Base.Math.asech ── メソッド
asech(A::AbstractMatrix)

行列ハイパボリックセカントの逆関数を正方行列 A に適用します。

Base.Math.acsch ── メソッド
acsch(A::AbstractMatrix)

行列ハイパボリックコセカントの逆関数を正方行列 A に適用します。

Base.Math.acoth ── メソッド
acoth(A::AbstractMatrix)

行列ハイパボリックコタンジェントの逆関数を正方行列 A に適用します。

LinearAlgebra.lyap ── 関数
lyap(A, C)

連続リアプノフ方程式 AX + XA' + C = 0 の解 X を計算します。A には実部がゼロの固有値が存在せず、どの二つの固有値も互いの共役の符号を反転させたものでないとします。

julia> A = [3. 4.; 5. 6]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> B = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> X = lyap(A, B)
2×2 Array{Float64,2}:
  0.5  -0.5
 -0.5   0.25

julia> A*X + X*A' + B
2×2 Array{Float64,2}:
 0.0          6.66134e-16
 6.66134e-16  8.88178e-16
LinearAlgebra.sylvester ── 関数
sylvester(A, B, C)

シルベスター方程式 AX + XB + C = 0 の解 X を計算します。A, B, C は計算が成立する次元を持ち、A-B には実部が等しい固有値が存在しないとします。

julia> A = [3. 4.; 5. 6]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> B = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> C = [1. 2.; -2. 1]
2×2 Array{Float64,2}:
  1.0  2.0
 -2.0  1.0

julia> X = sylvester(A, B, C)
2×2 Array{Float64,2}:
 -4.46667   1.93333
  3.73333  -1.8

julia> A*X + X*B + C
2×2 Array{Float64,2}:
  2.66454e-15  1.77636e-15
 -3.77476e-15  4.44089e-16
LinearAlgebra.issuccess ── 関数
issuccess(F::Factorization)

行列の分解が成功したかどうかを判定します。

julia> F = cholesky([1 0; 0 1]);

julia> LinearAlgebra.issuccess(F)
true

julia> F = lu([1 0; 0 0]; check = false);

julia> LinearAlgebra.issuccess(F)
false
issymmetric(A) -> Bool

行列が対称かどうかを判定します。

julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> issymmetric(a)
true

julia> b = [1 im; -im 1]
2×2 Array{Complex{Int64},2}:
 1+0im  0+1im
 0-1im  1+0im

julia> issymmetric(b)
false
LinearAlgebra.isposdef ── 関数
isposdef(A) -> Bool

A のコレスキー分解の計算を試みることで、行列が (エルミートかつ) 正定値かどうかを判定します。isposdef! も参照してください。

julia> A = [1 2; 2 50]
2×2 Array{Int64,2}:
 1   2
 2  50

julia> isposdef(A)
true
LinearAlgebra.isposdef! ── 関数
isposdef!(A) -> Bool

A のコレスキー分解の計算を試みることで、行列が (エルミートかつ) 正定値かどうかを判定します。A は処理中に上書きされます。isposdef も参照してください。

julia> A = [1. 2.; 2. 50.];

julia> isposdef!(A)
true

julia> A
2×2 Array{Float64,2}:
 1.0  2.0
 2.0  6.78233
LinearAlgebra.istril ── 関数
istril(A::AbstractMatrix, k::Integer = 0) -> Bool

A が第 k 優対角要素およびそれより下にのみ非ゼロ要素を持つ下三角行列かどうかを判定します。

julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> istril(a)
false

julia> istril(a, 1)
true

julia> b = [1 0; -im -1]
2×2 Array{Complex{Int64},2}:
 1+0im   0+0im
 0-1im  -1+0im

julia> istril(b)
true

julia> istril(b, -1)
false
LinearAlgebra.istriu ── 関数
istriu(A::AbstractMatrix, k::Integer = 0) -> Bool

A が第 k 優対角要素およびそれより上にのみ非ゼロ要素を持つ上三角行列かどうかを判定します。

julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> istriu(a)
false

julia> istriu(a, -1)
true

julia> b = [1 im; 0 -1]
2×2 Array{Complex{Int64},2}:
 1+0im   0+1im
 0+0im  -1+0im

julia> istriu(b)
true

julia> istriu(b, 1)
false
LinearAlgebra.isdiag ── 関数
isdiag(A) -> Bool

行列が対角行列かどうかを判定します。

julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> isdiag(a)
false

julia> b = [im 0; 0 -im]
2×2 Array{Complex{Int64},2}:
 0+1im  0+0im
 0+0im  0-1im

julia> isdiag(b)
true
ishermitian(A) -> Bool

行列がエルミート行列かどうかを判定します。

julia> a = [1 2; 2 -1]
2×2 Array{Int64,2}:
 1   2
 2  -1

julia> ishermitian(a)
true

julia> b = [1 im; -im 1]
2×2 Array{Complex{Int64},2}:
 1+0im  0+1im
 0-1im  1+0im

julia> ishermitian(b)
true
Base.transpose ── 関数
transpose(A)

行列の転置です。遅延評価が行われます。transpose の返り値を改変すると A も適切に変更されます。遅延された転置のラッパーを表す Transpose 型の Transpose(A) を多くの場合で返しますが、いつもそうだとは限りません。この操作は再帰的なことに注意してください。

この操作は線形代数の文脈での利用が意図されています。一般的なデータ操作では再帰的でない permutedims を使ってください。

julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> transpose(A)
2×2 Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 3+2im  8+7im
 9+2im  4+6im
LinearAlgebra.transpose! ── 関数
transpose!(dest,src)

配列 src を転置して、結果を事前にアロケートされた配列 dest に格納します。dest(size(src,2),size(src,1)) のサイズを持つべきです。インプレースの転置はサポートされず、srcdest のメモリ領域が重なっていると予期せぬ結果が生じます。

julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> B = zeros(Complex{Int64}, 2, 2)
2×2 Array{Complex{Int64},2}:
 0+0im  0+0im
 0+0im  0+0im

julia> transpose!(B, A);

julia> B
2×2 Array{Complex{Int64},2}:
 3+2im  8+7im
 9+2im  4+6im

julia> A
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im
Transpose

別の場所にある線形代数オブジェクトに対する転置ビューを表す遅延ラッパー型です。転置されるオブジェクトは通常 AbstractVector または AbstractMatrix ですが、Factorization であることもあります。通常 Transpose コンストラクタを直接使うべきではなく、代わりに transpose 関数を使うべきです。ビューを具現化するには copy を使ってください。

この型は線形代数の文脈での利用が意図されています。一般的なデータ操作では permutedims を使ってください。

julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> transpose(A)
2×2 Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 3+2im  8+7im
 9+2im  4+6im
Base.adjoint ── 関数
A'
adjoint(A)

随伴 (共役転置) です。遅延評価が行われます。adjoint は要素に対して再帰的に適用されることに注意してください。

adjoint は数値型に対して複素共役を返します。そのため実数に対する adjoint は恒等関数です。

この操作は線形代数の文脈での利用が意図されています。一般的なデータ操作では permutedims を使ってください。

julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> adjoint(A)
2×2 Adjoint{Complex{Int64},Array{Complex{Int64},2}}:
 3-2im  8-7im
 9-2im  4-6im

julia> x = [3, 4im]
2-element Array{Complex{Int64},1}:
 3 + 0im
 0 + 4im

julia> x'x
25 + 0im
LinearAlgebra.adjoint! ── 関数
adjoint!(dest,src)

配列 src の共役転置を取り、結果を事前にアロケートされた配列 dest に格納します。dest(size(src,2),size(src,1)) のサイズを持つべきです。インプレースの随伴はサポートされず、srcdest のメモリ領域が重なっていると予期せぬ結果が生じます。

julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> B = zeros(Complex{Int64}, 2, 2)
2×2 Array{Complex{Int64},2}:
 0+0im  0+0im
 0+0im  0+0im

julia> adjoint!(B, A);

julia> B
2×2 Array{Complex{Int64},2}:
 3-2im  8-7im
 9-2im  4-6im

julia> A
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im
Adjoint

別の場所にある線形代数オブジェクトに対する随伴ビューを表す遅延ラッパー型です。転置されるオブジェクトは通常 AbstractVector または AbstractMatrix ですが、Factorization であることもあります。通常 Adjoint コンストラクタを直接使うべきではなく、代わりに adjoint 関数を使うべきです。ビューを具現化するには copy を使ってください。

この型は線形代数の文脈での利用が意図されています。一般的なデータ操作では permutedims を使ってください。

julia> A = [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> adjoint(A)
2×2 Adjoint{Complex{Int64},Array{Complex{Int64},2}}:
 3-2im  8-7im
 9-2im  4-6im
Base.copy ── メソッド
copy(A::Transpose)
copy(A::Adjoint)

遅延された行列の転置/随伴を積極評価します。転置/随伴は要素に対して再帰的に適用されることに注意してください。

この操作は線形代数の文脈での利用が意図されています。一般的なデータ操作では再帰的でない permutedims を使ってください。

julia> A = [1 2im; -3im 4]
2×2 Array{Complex{Int64},2}:
 1+0im  0+2im
 0-3im  4+0im

julia> T = transpose(A)
2×2 Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0-3im
 0+2im  4+0im

julia> copy(T)
2×2 Array{Complex{Int64},2}:
 1+0im  0-3im
 0+2im  4+0im
LinearAlgebra.stride1 ── 関数
stride1(A) -> Int

一つ目の次元で隣り合う添え字を持つ配列の要素がどれだけ離れているかを返します。要素型のサイズが一単位となります。

julia> A = [1,2,3,4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> LinearAlgebra.stride1(A)
1

julia> B = view(A, 2:2:4)
2-element view(::Array{Int64,1}, 2:2:4) with eltype Int64:
 2
 4

julia> LinearAlgebra.stride1(B)
2
LinearAlgebra.checksquare(A)

行列が正方であることを確認し、二つの次元に共通のサイズを返します。複数の引数に対してはベクトルを返します。

julia> A = fill(1, (4,4)); B = fill(1, (5,5));

julia> LinearAlgebra.checksquare(A, B)
2-element Array{Int64,1}:
 4
 5
LinearAlgebra.peakflops ── 関数
LinearAlgebra.peakflops(n::Integer=2000; parallel::Bool=false)

倍精度の gemm! を使ってマシンのピーク浮動小数点数演算 (FLOP) レートを計算します。引数が指定されないと、n = 2000 に対する n×n 行列の乗算を行います。内部の BLAS が複数のスレッドを使用するなら、FLOP レートはその分高くなります。BLAS が使用するスレッド数は BLAS.set_num_threads(n) で設定できます。

キーワード引数 paralleltrue に設定されると、peakflops は全てのワーカープロセッサで並列に実行されます。並列に実行されるときは一つの BLAS スレッドだけが使われ、並列計算全体の FLOP レートが返ります。このとき n は各プロセッサで解かれる問題のサイズを指します。

Julia 1.1

この関数は Julia 1.1 以降のバージョンでサポートされます。Julia 1.0 では標準ライブラリの InteractiveUtils モジュールから利用できます。

低水準行列演算

多くの行列演算にはインプレースなバージョンが存在し、そういった関数には事前にアロケートされたベクトルあるいは行列を渡すことができます。これはループで反復されるアロケートのオーバーヘッドを避けてクリティカルなコードを最適化するときに有用です。Julia の慣習に従い、インプレースな演算は名前の最後に ! が付きます (mul! など)。

LinearAlgebra.mul! ── 関数
mul!(Y, A, B) -> Y

行列-行列積または行列-ベクトル積 \(AB\) を計算し、結果を Y に格納します。Y に元々あった値は上書きされます。YA または B の別名であってはいけません。

julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; Y = similar(B); mul!(Y, A, B);

julia> Y
2×2 Array{Float64,2}:
 3.0  3.0
 7.0  7.0

実装

独自の行列型とベクトル型では、可能な限り三引数の mul! ではなく五引数の mul! を直接実装することが推奨されます。

mul!(C, A, B, α, β) -> C

インプレースな行列-行列または行列-ベクトルの積和演算 \(A B α + C β\) です。C は返り値で上書きされます。CA または B の別名であってはいけません。

Julia 1.3

五引数の mul! は Julia 1.3 以降でサポートされます。

julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; C=[1.0 2.0; 3.0 4.0];

julia> mul!(C, A, B, 100.0, 10.0) === C
true

julia> C
2×2 Array{Float64,2}:
 310.0  320.0
 730.0  740.0
LinearAlgebra.lmul! ── 関数
lmul!(a::Number, B::AbstractArray)

配列 B をスカラー a でスケールします。演算はインプレースに行われ、B は上書きされます。スカラーを右から乗じるには rmul! を使ってください。このスケール演算は aB の要素を * で乗算するときの意味論に従います。特に NaN±Inf のような有限でない数値に対しても * と同様の処理が行われます。

Julia 1.1

Julia 1.1 より前のバージョンでは、B に含まれる NaN±Inf は特別扱いされていました。

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

julia> lmul!(2, B)
2×2 Array{Int64,2}:
 2  4
 6  8

julia> lmul!(0.0, [Inf])
1-element Array{Float64,1}:
 NaN
lmul!(A, B)

行列-行列積 \(AB\) を計算し、結果を B に上書きし、B を返します。ここで A は特殊な行列型である必要があります。例えば Diagonal, UpperTriangular, LowerTriangular や何らかの直交行列を表す型です。QR も参照してください。

julia> B = [0 1; 1 0];

julia> A = LinearAlgebra.UpperTriangular([1 2; 0 3]);

julia> LinearAlgebra.lmul!(A, B);

julia> B
2×2 Array{Int64,2}:
 2  1
 3  0

julia> B = [1.0 2.0; 3.0 4.0];

julia> F = qr([0 1; -1 0]);

julia> lmul!(F.Q, B)
2×2 Array{Float64,2}:
 3.0  4.0
 1.0  2.0
LinearAlgebra.rmul! ── 関数
rmul!(A::AbstractArray, b::Number)

配列 A をスカラー b でスケールし、結果を A に上書きします。スカラーを左から乗じるには lmul! を使ってください。このスケール演算は A の要素と b* で乗算するときの意味論に従います。特に NaN±Inf のような有限でない数値に対しても * と同様の処理が行われます。

Julia 1.1

Julia 1.1 より前のバージョンでは、B に含まれる NaN±Inf は特別扱いされていました。

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

julia> rmul!(A, 2)
2×2 Array{Int64,2}:
 2  4
 6  8

julia> rmul!([NaN], 0.0)
1-element Array{Float64,1}:
 NaN
rmul!(A, B)

行列-行列積 \(AB\) を計算し、結果を B に上書きし、B を返します。ここで A は特殊な行列型である必要があります。例えば Diagonal, UpperTriangular, LowerTriangular や何らかの直交行列を表す型です。QR も参照してください。

julia> A = [0 1; 1 0];

julia> B = LinearAlgebra.UpperTriangular([1 2; 0 3]);

julia> LinearAlgebra.rmul!(A, B);

julia> A
2×2 Array{Int64,2}:
 0  3
 1  2

julia> A = [1.0 2.0; 3.0 4.0];

julia> F = qr([0 1; -1 0]);

julia> rmul!(A, F.Q)
2×2 Array{Float64,2}:
 2.0  1.0
 4.0  3.0
LinearAlgebra.ldiv! ── 関数
ldiv!(Y, A, B) -> Y

A \ B をインプレースに計算し、結果を Y に格納します。計算結果 Y を返します。

引数 A には 行列を指定するべきではありませんA には行列ではなく分解オブジェクト (例えば factorizecholesky の返り値) を使ってください。行列の分解は時間のかかる処理であり、さらに通常メモリをアロケートするためです (ただし lu! などを使えばインプレースな分解も可能です)。ldiv! が必要になるパフォーマンスクリティカルな状況では、通常 A の分解についても細かな調整が必要になります。

julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];

julia> X = [1; 2.5; 3];

julia> Y = zero(X);

julia> ldiv!(Y, qr(A), X);

julia> Y
3-element Array{Float64,1}:
  0.7128099173553719
 -0.051652892561983674
  0.10020661157024757

julia> A\X
3-element Array{Float64,1}:
  0.7128099173553719
 -0.05165289256198333
  0.10020661157024785
ldiv!(A, B)

A \ B をインプレースに計算し、結果を B に上書きします。

引数 A には 行列を指定するべきではありませんA には行列ではなく分解オブジェクト (例えば factorizecholesky の返り値) を使ってください。行列の分解は時間のかかる処理であり、さらに通常メモリをアロケートするためです (ただし lu! などを使えばインプレースな分解も可能です)。ldiv! が必要になるパフォーマンスクリティカルな状況では、通常 A の分解についても細かな調整が必要になります。

julia> A = [1 2.2 4; 3.1 0.2 3; 4 1 2];

julia> X = [1; 2.5; 3];

julia> Y = copy(X);

julia> ldiv!(qr(A), X);

julia> X
3-element Array{Float64,1}:
  0.7128099173553719
 -0.051652892561983674
  0.10020661157024757

julia> A\Y
3-element Array{Float64,1}:
  0.7128099173553719
 -0.05165289256198333
  0.10020661157024785
ldiv!(a::Number, B::AbstractArray)

配列 B の各要素をスカラー a で割り、結果を B に上書きします。右からスカラー除算を行うには rdiv! を使ってください。

julia> B = [1.0 2.0; 3.0 4.0]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> ldiv!(2.0, B)
2×2 Array{Float64,2}:
 0.5  1.0
 1.5  2.0
LinearAlgebra.rdiv! ── 関数
rdiv!(A, B)

A / B をインプレースに計算し、結果を A に上書きします。

引数 B には 行列を指定するべきではありませんB には行列ではなく分解オブジェクト (例えば factorizecholesky の返り値) を使ってください。行列の分解は時間のかかる処理であり、さらに通常メモリをアロケートするためです (ただし lu! などを使えばインプレースな分解も可能です)。rdiv! が必要になるパフォーマンスクリティカルな状況では、通常 B の分解についても細かな調整が必要になります。

rdiv!(A::AbstractArray, b::Number)

配列 A の各要素をスカラー b で割り、結果を A に上書きします。左からスカラー除算を行うには ldiv! を使ってください。

julia> A = [1.0 2.0; 3.0 4.0]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> rdiv!(A, 2.0)
2×2 Array{Float64,2}:
 0.5  1.0
 1.5  2.0

BLAS 関数

Julia における密な線形代数演算は (多くの科学計算と同様) LAPACK ライブラリを利用します。そして LAPACK は BLAS と呼ばれる線形代数演算の基礎単位を使って作られています。高度に最適化された BLAS 実装は全てのコンピューターアーキテクチャで利用可能であり、高性能な線形代数ルーチンでは BLAS 関数を直接呼び出すのが望ましいこともあります。

LinearAlgebra.BLAS は一部の BLAS 関数に対するラッパーを提供します。入力の配列のいずれかを変更する BLAS 関数は名前が '!' で終わります。通常 BLAS 関数には四つのメソッドが定義され、それぞれ要素型が Float64, Float32, ComplexF64, ComplexF32 の配列に対応します。

BLAS の文字引数

多くの BLAS 関数は、引数を転置すべきかどうか (trans/tX)、 どちらの三角部分を参照するか (uplo/ul)、三角行列の対角要素が全て 1 とみなされるかどうか (dA)、入力引数は行列乗算のどちら側に置かれるか (side) を示す文字引数を受け取ります。各引数に指定できる値は次の通りです:

行列の転置

trans/tX 意味
'N' 入力の行列 X は転置も随伴 (共役転置) もされない。
'T' 入力の行列 X は転置される。
'C' 入力の行列 X は随伴 (共役転置) される。

三角部分の参照

uplo/ul 意味
'U' 行列の三角部分だけが使われる。
'L' 行列の三角部分だけが使われる。

単位対角要素

diag/dX 意味
'N' 行列 X の対角要素が読み込まれる。
'U' 行列 X の対角要素は全て 1 とみなされる。

乗算の順序

side 意味
'L' 引数は行列-行列積の側に置かれる。
'R' 引数は行列-行列積の側に置かれる。
LinearAlgebra.BLAS ── モジュール

BLAS サブルーチンへのインターフェースです。

LinearAlgebra.BLAS.dot ── 関数
dot(n, X, incx, Y, incy)

配列 X から歩長を incx として取った n 要素と配列 Y から歩長を incy として取った n 要素からなる二つのベクトルのドット積です。

julia> BLAS.dot(10, fill(1.0, 10), 1, fill(1.0, 20), 2)
10.0
LinearAlgebra.BLAS.dotu ── 関数
dotu(n, X, incx, Y, incy)

配列 X から歩長を incx として取った n 要素と配列 Y から歩長を incy として取った n 要素からなる二つの複素ベクトルのドット積です。

julia> BLAS.dotu(10, fill(1.0im, 10), 1, fill(1.0+im, 20), 2)
-10.0 + 10.0im
LinearAlgebra.BLAS.dotc ── 関数
dotc(n, X, incx, U, incy)

配列 X から歩長を incx として取った n 要素と配列 Y から歩長を incy として取った n 要素からなる二つの複素ベクトルのドット積です。一つ目のベクトルの共役が取られます。

julia> BLAS.dotc(10, fill(1.0im, 10), 1, fill(1.0+im, 20), 2)
10.0 - 10.0im
blascopy!(n, X, incx, Y, incy)

配列 X から歩長を incx として取った n 要素を配列 Y に歩長を incy としてコピーします。Y を返します。

LinearAlgebra.BLAS.nrm2 ── 関数
nrm2(n, X, incx)

配列 X から歩長を incx として取った n 要素からなるベクトルの 2 ノルムを返します。

julia> BLAS.nrm2(4, fill(1.0, 8), 2)
2.0

julia> BLAS.nrm2(1, fill(1.0, 8), 2)
1.0
LinearAlgebra.BLAS.asum ── 関数
asum(n, X, incx)

配列 X から歩長を incx として取った最初の n 要素の絶対値の和を返します。

julia> BLAS.asum(5, fill(1.0im, 10), 2)
5.0

julia> BLAS.asum(2, fill(1.0im, 10), 5)
2.0
LinearAlgebra.axpy! ── 関数
axpy!(a, X, Y)

YX*a + Y で上書きします。a はスカラーです。Y を返します。

julia> x = [1; 2; 3];

julia> y = [4; 5; 6];

julia> BLAS.axpy!(2, x, y)
3-element Array{Int64,1}:
  6
  9
 12
LinearAlgebra.axpby! ── 関数
axpby!(a, X, b, Y)

YX*a + Y*b で上書きします。ab はスカラーです。Y を返します。

julia> x = [1., 2, 3];

julia> y = [4., 5, 6];

julia> BLAS.axpby!(2., x, 3., y)
3-element Array{Float64,1}:
 14.0
 19.0
 24.0
LinearAlgebra.BLAS.scal! ── 関数
scal!(n, a, X, incx)

配列 X から歩長を incx として取った最初の n 要素を、a でスケールした値で上書きします。X を返します。

LinearAlgebra.BLAS.scal ── 関数
scal(n, a, X, incx)

配列 X から歩長を incx として取った最初の n 要素を a でスケールした配列を返します。

LinearAlgebra.BLAS.iamax ── 関数
iamax(n, dx, incx)
iamax(dx)

dx の絶対値が最大の要素の添え字を検索します。ndx の長さで、incxn の歩長です。nincx を指定しないと、n=length(dx)incx=stride1(dx) がデフォルトで使われます。

LinearAlgebra.BLAS.ger! ── 関数
ger!(alpha, x, y, A)

ベクトル xy による行列 A のランク 1 の更新 A = alpha*x*y' + A を行います。

LinearAlgebra.BLAS.syr! ── 関数
syr!(uplo, alpha, x, A)

対称行列 A に対してベクトル x によるランク 1 の更新 A = alpha*x*transpose(x) + A を行います。uploA のどちらの三角部分が更新されるかを制御します。A を返します。

LinearAlgebra.BLAS.syrk! ── 関数
syrk!(uplo, trans, alpha, A, beta, C)

対称行列 C に対するランク k の更新を行います。trans の値に応じて alpha*A*transpose(A) + beta*C または alpha*transpose(A)*A + beta*C を計算し、その結果で C を更新します。uplo が指定する C の三角部分だけが使われます。C を返します。

LinearAlgebra.BLAS.syrk ── 関数
syrk(uplo, trans, alpha, A)

trans の値に応じて alpha*A*transpose(A) または alpha*transpose(A)*A を計算し、計算結果の uplo が指定する三角部分を返します3

syr2k!(uplo, trans, alpha, A, B, beta, C)

対称行列 C に対してランク 2k の更新を行います。trans の値に応じて alpha*A*transpose(B) + alpha*B*transpose(A) + beta*C または alpha*transpose(A)*B + alpha*transpose(B)*A + beta*C を計算し、その結果で C を更新します。uplo が指定する C の三角部分だけが使われます。C を返します。

LinearAlgebra.BLAS.syr2k ── 関数
syr2k(uplo, trans, alpha, A, B)

alpha*A*transpose(B) + alpha*B*transpose(A) または alpha*transpose(A)*B + alpha*transpose(B)*A を計算し、計算結果の uplo が指定する三角部分を返します。どちらを計算するかは trans によって指定されます。

syr2k(uplo, trans, A, B)

A*transpose(B) + B*transpose(A) または transpose(A)*B + transpose(B)*A を計算し、計算結果の uplo が指定する三角部分を返します。どちらを計算するかは trans によって指定されます。

LinearAlgebra.BLAS.her! ── 関数
her!(uplo, alpha, x, A)

複素数配列専用のメソッドです。エルミート行列 A に対するベクトル x によるランク 1 の更新 A = alpha*x*x' + A を行います。uploA のどちらの三角部分が更新されるかを制御します。A を返します。

LinearAlgebra.BLAS.herk! ── 関数
herk!(uplo, trans, alpha, A, beta, C)

複素数配列専用のメソッドです。エルミート行列 C に対するランク k の更新を行います。trans の値に応じて alpha*A*A' + beta*C または alpha*A'*A + beta*C を計算し、その結果で C を更新します。Cuplo が指定する三角部分だけが更新されます。C を返します。

LinearAlgebra.BLAS.herk ── 関数
herk(uplo, trans, alpha, A)

複素数配列専用のメソッドです。trans の値に応じて alpha*A*A' または alpha*A'*A を計算し、計算結果の uplo が指定する三角部分を返します。

her2k!(uplo, trans, alpha, A, B, beta, C)

エルミート行列 C にランク 2k の更新を行います。trans の値に応じて alpha*A*B' + alpha*B*A' + beta*C または alpha*A'*B + alpha*B'*A + beta*CC が更新されます。beta はスカラーである必要があります。Cuplo が指定する三角部分だけが使われます。C を返します。

LinearAlgebra.BLAS.her2k ── 関数
her2k(uplo, trans, alpha, A, B)

trans の値に応じて alpha*A*B' + alpha*B*A' または alpha*A'*B + alpha*B'*A を計算し、計算結果の uplo が指定する三角部分を返します。

her2k(uplo, trans, A, B)

trans の値に応じて A*B' + B*A' または A'*B + B'*A を計算し、計算結果の uplo が指定する三角部分を返します。

LinearAlgebra.BLAS.gbmv! ── 関数
gbmv!(trans, m, kl, ku, alpha, A, x, beta, y)

trans の値に応じて、ベクトル yalpha*A*x + beta*y または alpha*A'*x + beta*y で更新します。行列 A は次元 m×size(A,2) の一般帯行列であり、対角要素より下に kl 段の要素を持ち、対角要素より上に ku 段の要素を持ちます。alphabeta はスカラーです。更新された y を返します。

LinearAlgebra.BLAS.gbmv ── 関数
gbmv(trans, m, kl, ku, alpha, A, x)

trans の値に応じて alpha*A*x または alpha*A'*x を返します。行列 A は次元 m×size(A,2) の一般帯行列であり、対角要素より下に kl 段の要素を持ち、対角要素より上に ku 段の要素を持ちます。alpha はスカラーです。

LinearAlgebra.BLAS.sbmv! ── 関数
sbmv!(uplo, k, alpha, A, x, beta, y)

ベクトル yalpha*A*x + beta*y で更新します。A はオーダー size(A,2) の対称帯行列であり、対角要素から k 段上までの優対角要素が引数 A に格納されます。A のストレージレイアウトは http://www.netlib.org/lapack/explore-html/ の Modules, Reference BLAS, Level2 BLAS で説明されています。Auplo が指定する三角部分だけが使われます。

更新された y を返します。

LinearAlgebra.BLAS.sbmv ── メソッド
sbmv(uplo, k, alpha, A, x)

alpha*A*x を返します。A はオーダー size(A,2) の対称帯行列で、対角要素から k 段上までの優対角要素が引数 A に格納されます。Auplo が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.sbmv ── メソッド
sbmv(uplo, k, A, x)

A*x を返します。A はオーダー size(A,2) の対称帯行列であり、対角要素から k 段上までの優対角要素が引数 A に格納されます。Auplo が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.gemm! ── 関数
gemm!(tA, tB, alpha, A, B, beta, C)

alpha*A*B + beta*C もしくは tAtB の値に応じて定まる他の三つの値のどれかで C を更新します。更新された C を返します。

LinearAlgebra.BLAS.gemm ── メソッド
gemm(tA, tB, alpha, A, B)

alpha*A*B もしくは tAtB の値に応じて定まる他の三つの値のどれかを返します。更新された C を返します。

LinearAlgebra.BLAS.gemm ── メソッド
gemm(tA, tB, A, B)

A*B もしくは tAtB の値に応じて定まる他の三つの値のどれかを返します。更新された C を返します。

LinearAlgebra.BLAS.gemv! ── 関数
gemv!(tA, alpha, A, x, beta, y)

tA の値に応じて alpha*A*x + beta*y または alpha*A'x + beta*y を計算し、その値で y を更新します。alphabeta はスカラーです。更新された y を返します。

LinearAlgebra.BLAS.gemv ── メソッド
gemv(tA, alpha, A, x)

tA の値に応じて alpha*A*x または alpha*A'x を返します。alphabeta はスカラーです。

LinearAlgebra.BLAS.gemv ── メソッド
gemv(tA, A, x)

tA の値に応じて A*x または A'x を返します。

LinearAlgebra.BLAS.symm! ── 関数
symm!(side, ul, alpha, A, B, beta, C)

side の値に応じて alpha*A*B + beta*C または alpha*B*A + beta*CC を更新します。A は対称とみなされ、ul が指定する三角部分だけが使われます。更新された C を返します。

LinearAlgebra.BLAS.symm ── メソッド
symm(side, ul, alpha, A, B)

side の値に応じて alpha*A*B または alpha*B*A を返します。A は対称とみなされ、ul が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.symm ── メソッド
symm(side, ul, A, B)

side の値に応じて A*B または B*A を返します。A は対称とみなされ、ul が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.symv! ── 関数
symv!(ul, alpha, A, x, beta, y)

ベクトル yalpha*A*x + beta*y で更新します。A は対称とみなされ、ul が指定する三角部分だけが使われます。alphabeta はスカラーです。更新された y を返します。

LinearAlgebra.BLAS.symv ── メソッド
symv(ul, alpha, A, x)

alpha*A*x を返します。A は対称とみなされ、ul が指定する三角部分だけが使われます。alpha はスカラーです。

LinearAlgebra.BLAS.symv ── メソッド
symv(ul, A, x)

A*x を返します。A は対称とみなされ、ul が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.hemm! ── 関数
hemm!(side, ul, alpha, A, B, beta, C)

side の値に応じて alpha*A*B + beta*C または alpha*B*A + beta*CC を更新します。A はエルミートとみなされ、ul が指定する三角部分だけが使われます。更新された C を返します。

LinearAlgebra.BLAS.hemm ── メソッド
hemm(side, ul, alpha, A, B)

side の値に応じて alpha*A*B または alpha*B*A を返します。A はエルミートとみなされ、ul が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.hemm ── メソッド
hemm(side, ul, A, B)

side の値に応じて A*B または B*A を計算します。A はエルミートとみなされ、ul が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.hemv! ── 関数
hemv!(ul, alpha, A, x, beta, y)

ベクトル yalpha*A*x + beta*y で更新します。A はエルミートとみなされ、ul が指定する三角部分だけが使われます。alphabeta はスカラーです。更新された y を返します。

LinearAlgebra.BLAS.hemv ── メソッド
hemv(ul, alpha, A, x)

alpha*A*x を返します。A はエルミートとみなされ、ul が指定する三角部分だけが使われます。alpha はスカラーです。

LinearAlgebra.BLAS.hemv ── メソッド
hemv(ul, A, x)

A*x を返します。A はエルミートとみなされ、ul が指定する三角部分だけが使われます。

LinearAlgebra.BLAS.trmm! ── 関数
trmm!(side, ul, tA, dA, alpha, A, B)

alpha*A*B もしくは sidetA の値で定まる他の三つの値のどれかで B を更新します。Aul が指定する三角部分だけが使われます。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。更新された B を返します。

LinearAlgebra.BLAS.trmm ── 関数
trmm(side, ul, tA, dA, alpha, A, B)

alpha*A*B もしくは sidetA の値で定まる他の三つの値のどれかを返します。Aul が指定する三角部分だけが使われます。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。

LinearAlgebra.BLAS.trsm! ── 関数
trsm!(side, ul, tA, dA, alpha, A, B)

A*X = alpha*B もしくは sidetA の値で定まる他の三つの方程式のどれかの解 XB に上書きします。Aul で定まる三角部分だけが使われます。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。更新された B を返します。

LinearAlgebra.BLAS.trsm ── 関数
trsm(side, ul, tA, dA, alpha, A, B)

A*X = alpha*B あるいは sidetA の値で定まる他の三つの方程式のどれかの解 X を返します。Aul で定まる三角部分だけが使われます。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。

LinearAlgebra.BLAS.trmv! ── 関数
trmv!(ul, tA, dA, A, b)

op(A)*b を返します。optA で定まります。Aul が指定する三角部分だけが使われます。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。乗算は b でインプレースに計算されます。

LinearAlgebra.BLAS.trmv ── 関数
trmv(ul, tA, dA, A, b)

op(A)*b を返します。optA で定まります。Aul が指定する三角部分だけが使われます。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。

LinearAlgebra.BLAS.trsv! ── 関数
trsv!(ul, tA, dA, A, b)

A*x = b もしくは tAul の値で定まる他の三つの方程式のどれかの解を b に上書きします。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。更新された b を返します。

LinearAlgebra.BLAS.trsv ── 関数
trsv(ul, tA, dA, A, b)

A*x = b もしくは tAul の値で定まる他の三つの方程式のどれかの解を返します。dAA の対角要素を読み込むか、それとも全て 1 とみなすかを定めます。

set_num_threads(n)

BLAS ライブラリが使用するスレッド数を設定します。

LAPACK 関数

LinearAlgebra.LAPACK は線形代数に関連する LAPACK 関数の一部に対するラッパーを提供します。このモジュールに含まれる関数で入力を上書きするものには、名前の最後に ! が付いています。

このモジュールの関数には通常四つのメソッドが定義されており、それぞれ要素型が Float64, Float32, ComplexF64, ComplexF32 の配列に対応します。

Julia が提供する LAPACK API は将来変更されることに注意してください。この API はユーザーが直接利用するためのものではないので、現在この API に含まれる関数の将来のリリースにおけるサポート/非推奨化に関しては何も保証されません。

LinearAlgebra.LAPACK ── モジュール

LAPACK サブルーチンに対するインターフェースです。

gbtrf!(kl, ku, m, AB) -> (AB, ipiv)

帯行列 AB の LU 分解を計算します。kl は対角要素より下にある非ゼロな帯の個数、ku は対角要素より上にある非ゼロな帯の個数、mAB の一つ目の次元のサイズを表します。インプレースに計算された LU 分解と使われたピボットのベクトル ipiv を返します。

gbtrs!(trans, kl, ku, m, AB, ipiv, B)

方程式 AB * X = B を解きます。transAB を置換するかどうかを定めます。trans に指定できる値は N (転置しない)・T (転置する)・C (共役転置する) のいずれかです。kl は対角要素より下にある非ゼロな帯の個数、ku は対角要素より上にある非ゼロな帯の個数、mAB の一つ目の次元のサイズ、ipivgbtrf! が返すピボットのベクトルを表します。計算結果のベクトルあるいは行列で B を上書きし、それを返します。

gebal!(job, A) -> (ilo, ihi, scale)

固有系やシューア分解の計算の前処理として、A のバランスを調整します。job に指定できる値は N (置換/スケールを行わない)・P (置換を行う)・S (スケールを行う)・B(置換/スケールを両方行う)のいずれかです。A をインプレースに改変して ilo, ihi, scale を返します。置換が有効だと、j > i かつ「1 < j < ilo または j > ihi」のとき A[i,j] = 0 が成り立ちます。scale には行われたスケーリング/置換に関する情報が格納されます。

gebak!(job, side, ilo, ihi, scale, V)

gebal! を使ってバランスを調整された行列の固有ベクトル V に逆方向のスケーリングと置換を施すことで、元の行列の固有ベクトルに変換します。side に指定できる値は L (左固有ベクトルが変形される) と R (右固有ベクトルが変換される) です。

gebrd!(A) -> (A, d, e, tauq, taup)

A をインプレースに二重対角形式 A = QBP' へ簡約します。返り値は二重対角行列 B に書き換わった A, B の対角要素を含む d, B の非対角要素を含む e, Q を表す基本鏡映子係数からなる tauq, P を表す基本鏡映子係数からなる taup です。

gelqf!(A, tau)

A の LQ 分解 A = LQ を計算します。tau には分解の基本鏡映子に対するパラメータを表すスカラーが格納されます。tauA の小さい方の次元以上の長さを持つ必要があります。

インプレースに改変した Atau を返します。

gelqf!(A) -> (A, tau)

A の LQ 分解 A = LQ を計算します。

A をインプレースに改変し、Atau を返します。tau は分解の基本鏡映子に対するパラメータを表すスカラーが格納されたベクトルです。

geqlf!(A, tau)

A の QL 分解 A = QL を計算します。tau には分解の基本鏡映子に対するパラメータを表すスカラーが格納されます。tauA の小さい方の次元以上の長さを持つ必要があります。

インプレースに改変した Atau を返します。

geqlf!(A) -> (A, tau)

A の QL 分解 A = QL を計算します。

A をインプレースに改変し、Atau を返します。tau は分解の基本鏡映子に対するパラメータを表すスカラーが格納されたベクトルです。

geqrf!(A, tau)

A の QR 分解 A = QR を計算します。tau には分解の基本鏡映子に対するパラメータを表すスカラーが格納されます。tauA の小さい方の次元以上の長さを持つ必要があります。

インプレースに改変した Atau を返します。

geqrf!(A) -> (A, tau)

AQR 分解 A = QR を計算します。

A をインプレースに改変し、Atau を返します。tau は分解の基本鏡映子に対するパラメータを表すスカラーが格納されたベクトルです。

geqp3!(A, [jpvt, tau]) -> (A, tau, jpvt)

レベル 3 BLAS を使ってピボット選択付き QR 分解 AP = QR を計算します。ピボット行列 Pjpvt で表され、tau には基本鏡映子が格納されます。jpvttau は事前にアロケートされた配列を渡すための省略可能な配列です。渡す場合には、jpvtA の二番目の次元以上の長さを持つ必要があり、tauA の小さい方の次元以上の長さを持つ必要があります。

A, jpvt, tau はインプレースに改変されます。

gerqf!(A, tau)

A の RQ 分解 A = RQ を計算します。tau には分解の基本鏡映子に対するパラメータを表すスカラーが格納されます。tauA の小さい方の次元以上の長さを持つ必要があります。

インプレースに改変した Atau を返します。

gerqf!(A) -> (A, tau)

A の RQ 分解 A = RQ を計算します。

A をインプレースに改変し、Atau を返します。tau は分解の基本鏡映子に対するパラメータを表すスカラーが格納されたベクトルです。

geqrt!(A, T)

A のブロック QR 分解 A = QR を計算します。T には分解の基本鏡映子に対するパラメータを表す上三角ブロック鏡映子が格納されます。T の一番目の次元のサイズがブロックサイズを決定し、この値は 1 から n (A の二番目の次元のサイズ) の間である必要があります。T の二番目の次元のサイズは A の小さい方の次元のサイズと一致する必要があります。

インプレースに改変された AT を返します。

geqrt!(A, nb) -> (A, T)

A のブロック QR 分解 A = QR を計算します。nb がブロックを設定し、この値は 1 から n (A の二番目の次元のサイズ) の間である必要があります。

A をインプレースに改変し、AT を返します。T には分解の基本鏡映子に対するパラメータを表す上三角ブロック鏡映子が格納されます。

geqrt3!(A, T)

A のブロック QR 分解 A = QR を再帰的に計算します。T には分解の基本鏡映子に対するパラメータを表す上三角ブロック鏡映子が格納されます。T の一番目の次元のサイズがブロックサイズを決定し、この値は 1 から n (A の二番目の次元のサイズ) の間である必要があります。T の二番目の次元のサイズは A の小さい方の次元のサイズと一致する必要があります。

インプレースに改変された AT を返します。

geqrt3!(A) -> (A, T)

A のブロック QR 分解 A = QR を再帰的に計算します。

A をインプレースに改変し、AT を返します。T には分解の基本鏡映子に対するパラメータを表す上三角ブロック鏡映子が格納されます。

getrf!(A) -> (A, ipiv, info)

ピボット選択付き LU 分解 PA = LU を計算します。

返り値はインプレースに改変した A, ピボットの情報 ipiv, 分解の状態を示す情報コード info です。分解が成功すれば info = 0 となり、U[i,i] が特異値であれば info = i となり、エラーがあれば info < 0 となります。

tzrzf!(A) -> (A, tau)

上台形行列 A を上三角の形へインプレースに変形します。返り値は Atau であり、tau は変形の基本鏡映子に対するパラメータを表すスカラーからなるベクトルです。

ormrz!(side, trans, A, tau, C)

tzrzf! によって生成される変形が表す Q を行列 C に乗じます。side の値に応じて乗算は左から (side = 'L', Q*C) もしくは右から (side = 'R', C*Q) となり、trans = 'C' だと Q は転置されます。

gels!(trans, A, B) -> (F, B, ssr)

線形方程式 A * X = B, transpose(A) * X = B, adjoint(A) * X = B のいずれかを QR 分解または LQ 分解を使って解きます。行列/ベクトル B はインプレースに改変されて解となります。A は自身の QR 分解または LQ 分解で上書きされます。trans に指定できる値は 'N' (改変なし)・'T' (転置)・'C' (共役転置) のいずれかです。gels! は最小ノルム/最小二乗解を求めます。A は優決定系でも劣決定系でも構いません。解は B を通して返されます。

gesv!(A, B) -> (B, A, ipiv)

正方行列 A に対する線形方程式 A * X = BA の LU 分解を使って解きます。A は自身の LU 分解で上書きされ、B は解 X で上書きされます。ipiv には A の LU 分解で使われるピボットに関する情報が格納されます。

getrs!(trans, A, ipiv, B)

正方行列 A に対する線形方程式 A * X = B, transpose(A) * X = B, adjoint(A) * X = B のいずれかを解きます。行列/ベクトル B は解でインプレースに改変されます。A には getrf! で計算された LU 分解、ipiv にはそのときのピボット情報を入力します。trans に指定できる値は 'N' (改変なし)・'T' (転置)・'C' (共役転置) です。

getri!(A, ipiv)

getrf! による LU 分解を利用して A の逆行列を計算します。ipivgetrf! が出力するピボット情報であり、A は同じ関数が出力する LU 分解です。AA が表す行列の逆行列で上書きされます。

gesvx!(fact, trans, A, AF, ipiv, equed, R, C, B)
       -> (X, equed, R, C, B, rcond, ferr, berr, work)

A の LU 分解を使って、線形方程式 A * X = B (trans = 'N')・transpose(A) * X = B (trans = 'T')・adjoint(A) * X = B (trans = 'C') のいずれかを解きます。fact に指定できる値は 'E', 'F', 'N' のいずれかです。fact'E' のとき A は平衡化 (equilibrate) されて AF にコピーされ、'F' のとき AFipiv が事前に計算された LU 分解からの入力となり、'N' のとき AAF にコピーされてから分解されます。fact = 'F' のとき equed には 'N', 'R', 'C', 'B' のいずれかを指定できます。'N'A が平衡化されていないことを表し、'R'A に左から Diagonal(R) を乗じたことを表し、'C'A に右から Diagonal(C) を乗じたことを表し、'B'A に左から Diagonal(R) および右から Diagonal(C) を乗じたことを表します。fact = 'F' かつ「equed = 'R' または equed = 'B'」なら、R の要素は全て正である必要があります。

返り値の X は線形方程式の解です。equedfactN でないときに使われる出力であり、行われた平衡化を表します。R は行平衡化の対角要素、C は列平衡化の対角要素です。trans = 'N' かつ「equed = 'R' または equed = 'B'」 のとき B は平衡化された形式 Diagonal(R)*B で上書きされ、「trans = 'T' または trans = 'C'」かつ「equed = 'C' または equed = 'B'」のとき BDiagonal(C)*B で上書きされます。rcond は平衡化後の逆条件数、ferrX に含まれる解ベクトルに対する前進エラーの上界、work はピボット増大係数の逆数です。

gesvx!(A, B)

平衡化と転置を行わない単純化されたバージョンの gesvx! です。

gelsd!(A, B, rcond) -> (B, rnk)

A * X = B の最小ノルム解を計算します。A の SVD 分解を求めて問題を分割統治する方法で計算されます。B は解 X で上書きされます。rcond より小さい特異値はゼロとみなされます。B に格納した解と A の実効ランク rnk を返します。

gelsy!(A, B, rcond) -> (B, rnk)

A * X = B の最小ノルム解を計算します。A の QR 分解を求めて問題を分割統治する方法で計算されます。B は解 X で上書きされます。rcond より小さい特異値はゼロとみなされます。B に格納した解と A の実効ランク rnk を返します。

gglse!(A, c, B, d) -> (X,res)

等式制約 B*x = d が付いた方程式 A*x = c を解きます。等式 ||c - A*x||^2 = 0 を使って計算されます。解 x と残差の二乗和が返ります。

geev!(jobvl, jobvr, A) -> (W, VL, VR)

A の固有系を求めます。jobvl = 'N' なら A の左固有ベクトルは計算されず、jobvr = 'N' なら右固有ベクトルは計算されません。jobvl = 'V' または jobvr = 'V' なら対応する固有ベクトルが計算されます。固有値 W, 右固有ベクトル VR, 左固有ベクトル VL を返します。

gesdd!(job, A) -> (U, S, VT)

分割統治のアプローチで A の固有値分解 A = U*S*V' を求めます。job = 'A' だと U の全ての列および V' の全ての行が計算されます。job = 'N' だと U の列と V' の行は計算されません。job = 'O' だと A が (薄い) U の列と (薄い) V' の行のどちらかで上書きされ、もう一方が U または VT として返ります4job = 'S' だと (薄い) U の列と (薄い) V' の行が個別に計算されて返されます。

gesvd!(jobu, jobvt, A) -> (U, S, VT)

A の特異値分解 A = U*S*V' を求めます。jobu = 'A' だと U の全ての列が計算され、jobvt = 'A' だと V' の全ての行が計算されます。jobu = 'N' だと U の列は計算されず、jobvt = 'N' だと V' の行は計算されません。jobu = O だと A が (薄い) U の列で上書きされ、jobvt = O だと A が (薄い) V' の行で上書きされます。jobu = S だと (薄い) U の列が個別に計算されて返され、jobvt = S だと (薄い) V' の行が個別に計算されて返されます。jobujobvt の両方を O とすることはできません。

U, S, Vt を返します。SA の特異値です。

ggsvd!(jobu, jobv, jobq, A, B) -> (U, V, Q, alpha, beta, k, l, R)

AB の一般化特異値分解 U'*A*Q = D1*R, V'*B*Q = D2*R を求めます。D1 は対角要素に alpha を持ち、D2 は対角要素に beta を持ちます。jobu = U だと直交/ユニタリ行列 U が計算され、jobv = 'V' だと直交/ユニタリ行列 V が計算され、jobq = Q だと直交/ユニタリ行列 Q が計算されます。jobu, jobv, jobqN だと、対応する行列は計算されません。この関数はバージョン 3.6.0 より前の LAPACK でのみ利用可能です。

ggsvd3!(jobu, jobv, jobq, A, B) -> (U, V, Q, alpha, beta, k, l, R)

AB の一般化特異値分解 U'*A*Q = D1*R, V'*B*Q = D2*R を求めます。D1 は対角要素に alpha を持ち、D2 は対角要素に beta を持ちます。jobu = U だと直交/ユニタリ行列 U が計算され、jobv = 'V' だと直交/ユニタリ行列 V が計算され、jobq = Q だと直交/ユニタリ行列 Q が計算されます。jobu, jobv, jobqN だと、対応する行列は計算されません。この関数はバージョン 3.6.0 以降の LAPACK でのみ利用可能です。

geevx!(balanc, jobvl, jobvr, sense, A)
       -> (A, w, VL, VR, ilo, ihi, scale, abnrm, rconde, rcondv)

行列に対するバランス調整を行いながら A の固有系を求めます。jobvl = 'N' だと A の左固有ベクトルは計算されず、jobvr = 'N' だと A の右固有ベクトルは計算されません。jobvl = 'V' あるいは jobvr = 'V' だと対応する固有ベクトルが計算されます。balanc = 'N' だとバランス調整は行われず、balanc = P だと A には置換だけが行われスケールされず、balanc = S だと A にはスケールだけが行われ置換はされず、balanc = B だと A は置換とスケールの両方が行われます。sense = 'N' だと逆条件数は計算されず、sense = E だと逆条件数は固有値に対してのみ計算され、sense = 'V' だと逆条件数は右固有ベクトルに対してのみ計算され、sense = B だと逆条件数は右固有ベクトルと固有値に対して計算されます。sense = E または sense = B のときは、右と左の固有ベクトルが計算される必要があります。

ggev!(jobvl, jobvr, A, B) -> (alpha, beta, vl, vr)

AB の一般化固有値分解を計算します。jobvl = 'N' だと左固有ベクトルは計算されず、jobvr = 'N' だと右固有ベクトルは計算されません。jobvl = 'V' または jobvl = 'V' だと、対応する固有ベクトルが計算されます。

gtsv!(dl, d, du, B)

方程式 A * X = B を解きます。A は三重対角行列で、dl が劣対角要素、d が対角要素、du が優対角要素を表します。

XB を上書きし、それを返します。

gttrf!(dl, d, du) -> (dl, d, du, du2, ipiv)

三重対角行列の LU 分解を求めます。dl が劣対角要素、d が対角要素、du が優対角要素を表します。

dl, d, du をインプレースに改変し、dl, d, du および第二優対角要素 du2 とピボットベクトル ipiv を返します。

gttrs!(trans, dl, d, du, du2, ipiv, B)

gttrf! で計算された LU 分解を使って方程式 A * X = B (trans = 'N')・transpose(A) * X = B (trans = 'T')・adjoint(A) * X = B (trans = 'C') のいずれかを解きます。B は解 X で上書きされます。

orglq!(A, tau, k = length(tau))

A に対する gelqf! の呼び出し結果を使って、LQ 分解の行列 Q を明示的に求めます。AQ で上書きされます。

orgqr!(A, tau, k = length(tau))

A に対する geqrf! の呼び出し結果を使って、QR 分解の行列 Q を明示的に求めます。AQ で上書きされます。

orgql!(A, tau, k = length(tau))

A に対する geqlf! の呼び出し結果を使って、QL 分解の行列 Q を明示的に求めます。AQ で上書きされます。

orgrq!(A, tau, k = length(tau))

A に対する gerqf! の呼び出し結果を使って、RQ 分解の行列 Q を明示的に求めます。AQ で上書きされます。

ormlq!(side, trans, A, tau, C)

side = 'L' のとき Q * C (trans = 'N')・transpose(Q) * C (trans = 'T')・adjoint(Q) * C (trans = 'C') のいずれかを計算し、side = 'R' のときこれを反転した乗算を計算します。ここで Qgelqf! で計算された A の LQ 分解の Q です。C は上書きされます。

ormqr!(side, trans, A, tau, C)

side = 'L' のとき Q * C (trans = 'N')・transpose(Q) * C (trans = 'T')・adjoint(Q) * C (trans = 'C') のいずれかを計算し、side = 'R' のときこれを反転した乗算を計算します。ここで Qgeqrf! で計算された A の QR 分解の Q です。C は上書きされます。

ormql!(side, trans, A, tau, C)

side = 'L' のとき Q * C (trans = 'N')・transpose(Q) * C (trans = 'T')・adjoint(Q) * C (trans = 'C') のいずれかを計算し、side = 'R' のときこれを反転した乗算を計算します。ここで Qgeqlf! で計算された A の QL 分解の Q です。C は上書きされます。

ormrq!(side, trans, A, tau, C)

side = 'L' のとき Q * C (trans = 'N')・transpose(Q) * C (trans = 'T')・adjoint(Q) * C (trans = 'C') のいずれかを計算し、side = 'R' のときこれを反転した乗算を計算します。ここで Qgerqf! で計算された A の RQ 分解の Q です。C は上書きされます。

gemqrt!(side, trans, V, T, C)

side = 'L' のとき Q * C (trans = 'N')・transpose(Q) * C (trans = 'T')・adjoint(Q) * C (trans = 'C') のいずれかを計算し、side = 'R' のときこれを反転した乗算を計算します。ここで Qgeqrt! で計算された A の RQ 分解の Q です。C は上書きされます。

posv!(uplo, A, B) -> (A, B)

方程式 A * X = B の解を求めます。A は対称またはエルミートな正定値行列です。uplo = 'U' だと A の上コレスキー分解が計算され、uplo = 'L' だと A の下コレスキー分解が計算されます。A は自身のコレスキー分解で上書きされ、B は解 X で上書きされます。

potrf!(uplo, A)

正定値行列 A コレスキー分解を計算します。uplo = 'U' なら上コレスキー分解が、uplo = 'L' なら下コレスキー分解が計算されます。A は上書きされ、A と情報コードが返ります。

potri!(uplo, A)

potrf! でコレスキー分解された正定値行列 A の逆行列を計算します。uplo = 'U' なら上コレスキー分解、uplo = 'L' なら下コレスキー分解が行われたものとして処理されます。

自身の逆行列で上書きされた A が返ります。

potrs!(uplo, A, B)

方程式 A * X = B の解を求めます。Apotrf! でコレスキー分解された対称/エルミートな正定値行列です。uplo = 'U' なら上コレスキー分解、uplo = 'L' なら下コレスキー分解が行われたものとして処理されます。B は解 X で上書きされます。

pstrf!(uplo, A, tol) -> (A, piv, rank, info)

正定値行列 A のピボット選択付きコレスキー分解を計算します。uplo = 'U' なら上コレスキー分解、uplo = 'L' なら下コレスキー分解が計算されます。tol はユーザーが設定できる許容値です。A は自身のコレスキー分解で上書きされます。

行列 A, ピボット piv, A の階数 rank, 情報コード info を返します。分解が成功すると info = 0 となり、A が不定値または階数不足だと info = i > 0 となります。

ptsv!(D, E, B)

対称/エルミート5正定値三重対角行列 A に対する方程式 A * X = B を解きます。DA の対角要素、EA の副対角要素です。B は解 X で上書きされ、B が返ります。

pttrf!(D, E)

対角要素 D と副対角要素 E を持つ対称/エルミート正定値対称三重対角行列の LDLt 分解を計算します。DE が上書きされて返ります。

pttrs!(D, E, B)

対角要素 D と副対角要素 E を持つ対称/エルミート正定値対称三重対角行列 A に対して、方程式 A * X = Bpttrf! による LDLt 分解を使って解きます。B は解 X で上書きされます。

trtri!(uplo, diag, A)

三角行列 A の逆行列を求めます。Auplo = 'U' なら上三角であり、uplo = 'L' なら下三角です。diag = 'N' だと A1 とは限らない対角要素を持つとみなされ、diag = 'U' だと A の対角要素は全て 1 とみなされます。A は自身の逆行列で上書きされます。

trtrs!(uplo, trans, diag, A, B)

三角行列 A に対して A * X = B (trans = 'N')・transpose(A) * X = B (trans = 'T')・adjoint(A) * X = B (trans = 'C') のいずれかを解きます。Auplo = 'U' なら上三角であり、uplo = 'L' なら下三角です。diag = 'N' だと A1 とは限らない対角要素を持つとみなされ、diag = 'U' だと A の対角要素は全て 1 とみなされます。B は解 X で上書きされます。

trcon!(norm, uplo, diag, A)

三角行列 A の逆条件数を求めます。Auplo = 'U' なら上三角であり、uplo = 'L' なら下三角です。diag = 'N' だと A1 とは限らない対角要素を持つとみなされ、diag = 'U' だと A の対角要素は全て 1 とみなされます。norm = 'I' だと条件数は無限ノルムで計算され、norm = 'O' あるいは norm = '1' だと条件数は 1 ノルムで計算されます。

trevc!(side, howmny, select, T, VL = similar(T), VR = similar(T))

上三角行列 T の固有系を求めます。side = 'R' だと右固有ベクトルが計算され、side = 'L' だと左固有ベクトルが計算され、side = 'B' だと両方が計算されます。howmny = 'A' だと全ての固有ベクトルが返り、howmny = 'B' だと全ての固有ベクトルを VLVR を使って逆変換したものが返り、howmny = 'S' だと select に対応する固有ベクトルが返ります。

trrfs!(uplo, trans, diag, A, B, X, Ferr, Berr) -> (Ferr, Berr)

side = 'L' なら、A * X = B (trans = 'N')・transpose(A) * X = B (trans = 'T')・adjoint(A) * X = B (trans = 'C') のいずれかの解に含まれる誤差を概算します。side = 'R' だと方程式に含まれる乗算が反転して X * A などとなります。この関数の入力は trtrs! の出力です。A は下三角行列であり、diag = 'N' なら A の対角要素は 1 でないとみなされ、diag = 'U' なら A の対角要素は 1 とみなされます。Ferr は前進誤差、Berr は後退誤差であり、どちらも要素ごとに計算されます。

stev!(job, dv, ev) -> (dv, Zmat)

dv を対角要素、ev を副対角要素に持つ対称三重対角行列の固有系を計算します。job = 'N' だと固有値だけが計算されて dv として返され、job = 'V' だと固有ベクトルも計算されて Zmat として返されます。

stebz!(range, order, vl, vu, il, iu, abstol, dv, ev) -> (dv, iblock, isplit)

dv を対角要素、ev を副対角要素に持つ対称三重対角行列の固有系を計算します。range = 'A' だと全ての固有値が計算され、range = 'V' だと半開区間 (vl, vu] に含まれる固有値だけが計算され、range = 'I' だと il 番目から iu 番目の固有値が計算されます。order = B だと固有値は一つのブロック内で並べ替えられ、order = E だと全てのブロックを通して並べ替えられます。abstol は収束判定で使う許容値を設定します。

stegr!(jobz, range, dv, ev, vl, vu, il, iu) -> (w, Z)

dv を対角要素、ev を副対角要素に持つ対称三重対角行列に対して、固有値 (jobz = 'N' のとき) または固有値と固有ベクトル (jobz = 'V' のとき) を計算します。range = 'A' だと全ての固有値が計算され、range = 'V' だと半開区間 (vl, vu] に含まれる固有値だけが計算され、range = 'I' だと il 番目から iu 番目の固有値が計算されます。固有値は w として返され、固有ベクトルは Z として返されます。

stein!(dv, ev_in, w_in, iblock_in, isplit_in)

dv を対角要素、ev_in を副対角要素に持つ対称三重対角行列の固有ベクトルを計算します。w_in は求める固有ベクトルに対応する固有値を指定し、iblock_inw_in に含まれる固有値に対応する部分行列を指定し、isplit_in は行列ブロックの分け方を指定します。

syconv!(uplo, A, ipiv) -> (A, work)

三角行列に分解された対称行列 A を二つの行列 LD に変換します。uplo = 'U' のとき A は上三角行列で、uplo = 'L' のとき A は下三角行列です。ipiv は三角行列への分解におけるピボットベクトルです。ALD で上書きされます。

sysv!(uplo, A, B) -> (B, A, ipiv)

対称行列 A に対する方程式 A * X = B の解を求めます。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。B は解 X で上書きされ、A はバンチ-カウフマン分解で上書きされます。ipiv は分解で使われたピボットの情報です。

sytrf!(uplo, A) -> (A, ipiv, info)

対称行列 A のバンチ-カウフマン分解を計算します。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。

分解で上書きされた A, ピボットベクトル ipiv, エラー情報を表す非負整数 info が返ります。info が正のとき A は特異であり、分解における info 番目の対角要素がちょうどゼロに等しくなります。

sytri!(uplo, A, ipiv)

sytrf! の結果を使って対称行列 A の逆行列を計算します。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。A は逆行列で上書きされます。

sytrs!(uplo, A, ipiv, B)

sytrf! の結果を使って対称行列 A に対する方程式 A * X = B を解きます。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。

hesv!(uplo, A, B) -> (B, A, ipiv)

エルミート行列 A に対する方程式 A * X = B の解を求めます。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。B は解 X で上書きされ、A は自身のバンチ-カウフマン分解で上書きされます。ipiv には分解に関するピボット情報が含まれます。

hetrf!(uplo, A) -> (A, ipiv, info)

エルミート行列 A のバンチ-カウフマン分解を計算します。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。

返り値は分解で上書きされた A, ピボットベクトル ipiv, エラーコードを表す非負整数 info です。info が正のとき A は特異であり、分解における info 番目の対角要素がちょうどゼロに等しくなります。

hetri!(uplo, A, ipiv)

sytrf! の結果を使ってエルミート行列 A の逆行列を計算します。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。A は逆行列で上書きされます。

hetrs!(uplo, A, ipiv, B)

sytrf! の結果を使ってエルミート行列 A に対する方程式 A * X = B を解きます。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。

syev!(jobz, uplo, A)

対称行列 A の固有値 (jobz = 'N' のとき) または固有値と固有ベクトル (jobz = 'V' のとき) を求めます。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。

syevr!(jobz, range, uplo, A, vl, vu, il, iu, abstol) -> (W, Z)

対称行列 A の固有値 (jobz = 'N' のとき) または固有値と固有ベクトル (jobz = 'V' のとき) を求めます。uplo = 'U' なら A には上三角部分が格納され、uplo = 'L' なら A には下三角部分が格納されます。range = 'A' だと全ての固有値が計算され、range = 'V' だと半開区間 (vl, vu] に含まれる固有値だけが計算され、range = 'I' だと il 番目から iu 番目の固有値が計算されます。abstol を設定すれば収束判定の許容値として使われます。

固有値は W として返され、固有ベクトルは Z として返されます。

sygvd!(itype, jobz, uplo, A, B) -> (w, A, B)

対称行列 A と対称正定値行列 B の一般化固有値 (jobz = 'N' のとき) または一般化固有値と一般化固有ベクトル (jobz = 'V' のとき) を求めます。uplo = 'U' なら AB の上三角部分が使われ、uplo = 'L' なら AB の下三角部分が使われます。itype = 1 なら A*x = lambda*B*x が解かれ、itype = 2 なら A*B*x = lambda*x が解かれ、itype = 3 なら B*A*x = lambda*x が解かれます。

bdsqr!(uplo, d, e_, Vt, U, C) -> (d, Vt, U, C)

d を対角要素、e_ を副対角要素に持つ二重対角行列の特異値分解を計算します。uplo = 'U' なら e_ は優対角要素であり、uplo = 'L' なら e_ は劣対角要素です。Q' * C を計算することもできます。

特異値は d として返り、行列 CQ' * C で上書きされます。

bdsdc!(uplo, compq, d, e_) -> (d, e, u, vt, q, iq)

d を対角要素、e_ を副対角要素に持つ二重対角行列の特異値分解を分割統治法で計算します。uplo = 'U' なら e_ は優対角要素であり、uplo = 'L' なら e_ は劣対角要素です。compq = 'N' だと特異値だけが計算され、compq = 'I' だと特異値と特異ベクトルが計算され、compq = 'P' だと特異値と特異ベクトルがコンパクト形式で計算されます。この関数は実数型に対してのみ動作します。

特異値は d として返ります。compq = 'P' ならコンパクトな特異ベクトルが iq として返ります。

gecon!(normtype, A, anorm)

行列 A の逆条件数を求めます。normtype = 'I' だと条件数は無限ノルムで計算され、normtype = 'O' あるいは normtype = '1' だと条件数は 1 ノルムで計算されます。Agetrf! の返り値、anormnormtype が指定するタイプの A のノルムである必要があります。

gehrd!(ilo, ihi, A) -> (A, tau)

行列 A をヘッセンベルグ形式に変換します。もし Agebal! によってバランスを調整されているなら、iloihigebal! の出力である必要があります。そうでなければ ilo = 1 かつ ihi = size(A,2) が必要です。tau には分解の基本鏡映子が含まれます。

orghr!(ilo, ihi, A, tau)

gehrd! の返り値から直交/ユニタリ行列 Q を明示的に求めます。ilo, ihi, A, taugehrd! の出力および入力に対応しなければなりません。

gees!(jobvs, A) -> (A, vs, w)

行列 A の固有値 (jobvs = 'N' のとき) または固有値とシューアベクトル (jobvs = 'V' のとき) を計算します。A はシューア形式で上書きされます。

返り値は上書きされた A, シューアベクトル vs, 固有値を含む w です。

gges!(jobvsl, jobvsr, A, B) -> (A, B, alpha, beta, vsl, vsr)

AB の一般化固有値と一般化シューア形式、そして左シューアベクトル (jobsvl = 'V' のとき) または右シューアベクトル (jobvsr = 'V' のとき) を返します。

一般化固有値は alphabeta に返され、左シューアベクトルは vsl に返され、右シューアベクトルは vsr に返されます。

trexc!(compq, ifst, ilst, T, Q) -> (T, Q)

行列のシューア分解を並べ替えます。compq = 'V' だとシューアベクトル Q が並べ替えられ、compq = 'N' だと改変されません。ifstilst がベクトルの並び替えを指定します。

trsen!(compq, job, select, T, Q) -> (T, Q, w, s, sep)

行列のシューア分解を並び替えを行い、指定された場合には逆条件数を計算します。job = 'N' だと逆条件数は計算されず、job = 'E' だと指定した固有値の集合に対する条件数だけが計算され、job = 'V' だと不変部分空間に対する条件数だけが計算され、job = 'B' だと指定した固有値の集合および部分空間に対する条件数が計算されます。compq = 'V' だとシューアベクトル Q が更新され、compq = 'N' だとシューアベクトルは変更されません。select が固有値を指定します。

返り値は T, Q, 並び替えられた固有値 w, 固有値の集合に対する条件数 s, 不変部分空間の条件数 sep です。

tgsen!(select, S, T, Q, Z) -> (S, T, alpha, beta, Q, Z)

一般化シューア分解のベクトルを並び替えます。select が各集合に含まれる固有値を指定します。

trsyl!(transa, transb, A, B, C, isgn=1) -> (C, scale)

シルベスター行列方程式 A * X +/- X * B = scale*C を解きます。ここで AB はどちらも擬上三角行列です。transa = 'N' だと A は変更されず、transa = 'T' だと A は転置され、transa = 'C' だと A は共役転置されます。transbB についても同様です。isgn = 1 だと方程式 A * X + X * B = scale * C が解かれ、isgn = n だと方程式 A * X - X * B = scale * C が解かれます。

Xscale を返します。C は解 X で上書きされます。

  • [Bischof1987] C Bischof and C Van Loan, "The WY representation for products of Householder matrices", SIAM J Sci Stat Comput 8 (1987), s2-s13. doi:10.1137/0908009
  • [Schreiber1989] R Schreiber and C Van Loan, "A storage-efficient WY representation for products of Householder transformations", SIAM J Sci Stat Comput 10 (1989), 53-57. doi:10.1137/0910005
  • [Bunch1977] J R Bunch and L Kaufman, Some stable methods for calculating inertia and solving symmetric linear systems, Mathematics of Computation 31:137 (1977), 163-179. http://www.ams.org/journals/mcom/1977-31-137/S0025-5718-1977-0428694-0/.
  • [issue8859] Issue 8859, "Fix least squares", https://github.com/JuliaLang/julia/pull/8859
  • [B96] Åke Björck, "Numerical Methods for Least Squares Problems", SIAM Press, Philadelphia, 1996, "Other Titles in Applied Mathematics", Vol. 51. doi:10.1137/1.9781611971484
  • [S84] G. W. Stewart, "Rank Degeneracy", SIAM Journal on Scientific and Statistical Computing, 5(2), 1984, 403-413. doi:10.1137/0905030
  • [KY88] Konstantinos Konstantinides and Kung Yao, "Statistical analysis of effective singular values in matrix rank determination", IEEE Transactions on Acoustics, Speech and Signal Processing, 36(5), 1988, 757-763. doi:10.1109/29.1585
  • [H05] Nicholas J. Higham, "The squaring and scaling method for the matrix exponential revisited", SIAM Journal on Matrix Analysis and Applications, 26(4), 2005, 1179-1193. doi:10.1137/090768539
  • [AH12] Awad H. Al-Mohy and Nicholas J. Higham, "Improved inverse scaling and squaring algorithms for the matrix logarithm", SIAM Journal on Scientific Computing, 34(4), 2012, C153-C169. doi:10.1137/110852553
  • [AHR13] Awad H. Al-Mohy, Nicholas J. Higham and Samuel D. Relton, "Computing the Fréchet derivative of the matrix logarithm and estimating the condition number", SIAM Journal on Scientific Computing, 35(4), 2013, C394-C410. doi:10.1137/120885991
  • [BH83] Åke Björck and Sven Hammarling, "A Schur method for the square root of a matrix", Linear Algebra and its Applications, 52-53, 1983, 127-140. doi:10.1016/0024-3795(83)80010-X
  • [AH16] Mary Aprahamian and Nicholas J. Higham, "Matrix Inverse Trigonometric and Inverse Hyperbolic Functions: Theory and Algorithms", MIMS EPrint: 2016.4. https://doi.org/10.1137/16M1057577

  1. 訳注: ここでは引数の値に応じて複数のアルゴリズムが切り替わるという意味で多重アルゴリズム (polyalgorithm) という言葉が使われている。Julia では多重ディスパッチにより引数の型に応じて呼ばれるメソッドが切り替わるが、この \ ではさらに引数の値に応じてアルゴリズムが (根本的に) 切り替わるということ。[return]

  2. 訳注: 英語版でこの次にあるヘッセンベルグ分解の説明は hessenberg 関数の解説と重複していたので削除した。[return]

  3. 訳注: 英語版の記述を実際の動作に合うよう修正した。[return]

  4. 訳注: 英語版の記述を実際の動作に合うよう修正した。[return]

  5. 訳注: 英語版の記述を実際の動作に合うよう修正した。以降いくつかの関数でも同様とした。[return]