◆ dgels()

優決定または劣決定系連立一次方程式 Ax = b の解 (フルランク)

目的

dgelsはm×n実行列Aあるいはその転置からなる優決定または劣決定系連立1次方程式をAのQRあるいはLQ分解を用いて解く. Aはフルランクであること.

下記オプションが提供される:

trans = 'N' かつ m >= n の場合: 優決定系の最小二乗解を求める. すなわち, 次の最小二乗問題を解く.
|| B - A*X || を最小化する.
trans = 'N' かつ m < n の場合: 劣決定系 A * X = B の最小ノルム解を求める.
trans = 'T' かつ m >= n の場合: 劣決定系 A^T * X = B の最小ノルム解を求める.
trans = 'T' かつ m < n の場合: 優決定系の最小二乗解を求める. すなわち, 次の最小二乗問題を解く.
|| B - A^T*X || を最小化する.

いくつかの右辺ベクトル b および解ベクトル x を1回の呼び出しで扱うことができる. これらのベクトルは, m×nrhs右辺行列Bおよびn×nrhs解行列Xの列として格納される.

戻り値: info (Int32)
= 0: 正常終了
= -1: 入力パラメータ trans の誤り (trans != 'T'および'N')
= -2: 入力パラメータ m の誤り (m < 0)
= -3: 入力パラメータ n の誤り (n < 0)
= -4: 入力パラメータ a の誤り
= -5: 入力パラメータ b の誤り
= -6: 入力パラメータ nrhs の誤り (nrhs < 0)
= i > 0: Aの分解結果の三角行列のi番目の対角要素が0になった. 従って, Aはフルランクではない. 最小二乗解は求めることができなかった.

引数

[in]	trans	= 'N': 連立一次方程式はAからなる. = 'T': 連立一次方程式はA^Tからなる.
[in]	m	行列 A の行数. (m >= 0) (m = 0 の場合, bにゼロベクトルを返す)
[in]	n	行列 A の列数. (n >= 0) (n = 0 の場合, bにゼロベクトルを返す)
[in,out]	a	2次元配列 (Float64, m x n) [in] m×n行列 A. [out] m >= n: aはdgeqrfが返すQR分解結果により上書きされる. m < n: aはdgelqfが返すLQ分解結果により上書きされる.
[in,out]	b	1 または 2次元配列 (Float64, max(m, n) または max(m, n) x nrhs) [in] 右辺ベクトルが列ごとに格納された行列 B. Bは, m×nrhs行列(trans = 'N'の場合), または, n×nrhs行列(trans = 'T'の場合)である. [out] info = 0の場合, bは列ごとに格納された解ベクトルにより上書きされる: trans = 'N'かつm >= n: bの行0〜n-1に最小二乗解ベクトルが入る. 各列の解の残差二乗和はその列の行n〜m-1の要素の二乗和で与えられる. trans = 'N'かつm < n: bの行0〜n-1に最小ノルム解ベクトルが入る. trans = 'T'かつm >= n: bの行0〜m-1に最小ノルム解ベクトルが入る. trans = 'T'かつm < n: bの行0〜m-1に最小二乗解ベクトルが入る. 各列の解の残差二乗和はその列の行m〜n-1の要素の二乗和で与えられる.
[in]	nrhs	(省略可) 右辺の数, すなわち, 行列BおよびXの列数. (nrhs >= 0) (nrhs = 0 の場合, 処理を行わずに戻る) (省略時 = 1)

使用例: 優決定系連立1次方程式 Ax = B の最小二乗解を求める. また, 分散を求める. ただし,
( -1.06 0.48 -0.04 )

A = ( -1.19 0.73 -0.24 )

( 1.97 -0.89 0.56 )

( 0.68 -0.53 0.08 )

( 0.3884 )

B = ( 0.1120 )

( -0.3644 )

( -0.0002 )

とする.
function TestDgels()

m = 4

n = 3

a = [ -1.06 0.48 -0.04;

-1.19 0.73 -0.24;

1.97 -0.89 0.56;

0.68 -0.53 0.08 ]

b = [ 0.3884, 0.1120, -0.3644, -0.0002 ]

info = dgels('N', m, n, a, b)

println(b[1], " ", b[2], " ", b[3], " ", info)

if info == 0

ci = Vector{Cdouble}(undef, n)

info = dgecov(0, n, a, ci)

println(ci, " ", info)

end

end

dgels
function dgels(trans::Char, m::Integer, n::Integer, a::Array{Float64}, b::Array{Float64}, nrhs::Integer=1)
優決定または劣決定系連立一次方程式 Ax = b の解 (フルランク)

dgecov
function dgecov(job::Integer, n::Integer, a::Array{Float64}, ci::Array{Float64})
線形最小二乗問題の分散・共分散行列 (dgels用)

実行結果: > TestDgels()

-0.8200000000000006 -0.9400000000000007 0.7400000000000008 0

[6.46959967542666, 16.735040821891918, 18.717753277322856] 0