---
title: Proc FCMPで行列計算
tags:  #[第11回大阪sas勉強会]  
author: [森岡裕[SASユーザー総会世話人]](https://convert.docswell.com/user/6484025)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/G75RWMKLE4.jpg?width=480
description: [第11回大阪sas勉強会]
published: April 15, 25
canonical: https://convert.docswell.com/s/6484025/K6VDQE-2025-04-15-043654
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/G75RWMKLE4.jpg)

Proc FCMPで行列計算
イーピーエス株式会社
森岡裕


# Page. 2

![Page Image](https://bcdn.docswell.com/page/9J2DQ9W3JR.jpg)

SAS® IMLはパワフルで柔軟な行列プログラミング言語であり，シンプルな構文と広
大な関数ライブラリで、容易に数式を統計プログラミング文に変換できます
原則として，行列計算は，IMLをライセンス導入していただいての実施が推奨です．
ただし，仮にIMLライセンスを持っていない場合でも行列計算を実施するための仕組
みがSASに備わっています．
具体的にはProc FCMPの行列操作ルーチンを利用する方法と，Proc DS2プロシジャ
のMATRIXパッケージを利用する方法です．
両方法とも，行列の積や逆行列，行列の転置などに必要な機能が準備されています
が，本発表ではFCMPの方を取り上げます．
理由としては，単位行列の生成，行列式の導出，コレズキー分解などの機能が備え
付けで用意されていることと，DS2プロシジャの場合は．基本的な文法がSASではな
くDS2言語になるため，少し学習コストが高いと思われるところによります


# Page. 3

![Page Image](https://bcdn.docswell.com/page/DEY2W4L8EM.jpg)

使える部品は以下．IMLに比べれば，機能の貧弱さは否めないが，最低限は揃っている
行列操作ルーチン
説明
関数
説明
CALL MULT
2つの入力行列の乗法積を計算
READ_ARRAY
CALL INV
逆行列を計算
データセットを
FCMP配列(行列)
に変換
CALL TRANSPOSE
行列の転置
WRITE_ARRAY
CALL IDENTITY
入力行列を単位行列に変換
FCMP配列(行列)
をデータセットに
変換
CALL DET
行列式を計算
CALL COL
コレスキー分解を計算
CALL ELEMMULT
elementwise乗算(アダマール積,要素ごとの積)を計算
CALL ADDMATRIX
elementwise乗加算(要素ごとの加算)を計算
CALL SUBTRACTMATRIX
element-wide減算(要素ごとの減算)を計算
CALL POWER
正方行列を指定されているスカラー値に累乗
CALL FILLMATRIX
すべての要素値を指定値と置き換え
CALL ZEROMATRIX
すべての要素値を0と置き換え
CALL EXPMATRIX
行列etAを 計算(入力行列A、乗数tの場合)
基本的にはFCMP内での配列を
行列とみなして，行列操作ルーチン
で操作して行列計算するイメージ．
データセットを読み込んで行列にす
ることも，その逆も1関数の簡単な
定義のみでできるため，使い勝手は
よい．
FCMPはユーザー定義の関数やサブ
ルーチンを作成するものなので，足
りたい部品は自作せよとのことなの
かなとか


# Page. 4

![Page Image](https://bcdn.docswell.com/page/VJN89Y4978.jpg)

行列操作ルーチン
説明
CALL MULT
2つの入力行列の乗法積を計算
成立要件：最初の入力行列の列数が2番目の行列の行数と同じである必要があります
proc fcmp;
array mat1[2,3] (1,2,4,5,8,9);
array mat2[3,2] (3,6,2,1,0,10);
第1引数にX, 第2引数にY, 第3引数にZ(結果格納行列)
array result[2,2];
call mult(mat1, mat2, result);
put result=;
quit;
1 2 4
mat1 =
5 8 9
3 6
mat2 = 2 1
0 10
行×列 ＝ 積
最初の行列の要素を横に，2個めの行列の要素を縦に掛け算して
足していけばいい
1×3+2×2 + 4×0 1×6+2×1 + 4×10
mat1 × mat2 =
=
5×3+8×2 + 9×0 5×6+8×1 + 9×10
7
48
31 128


# Page. 5

![Page Image](https://bcdn.docswell.com/page/YE9K2PQ3E3.jpg)

関数
説明
WRITE_ARRAY
FCMP配列(行列)をデータセットに変換
第一引数に”出力データセット名”, 第二引数に配列名，第三以降は変数名を指定の場合使用
proc fcmp;
array mat1[2,3] (1,2,4,5,8,9);
array mat2[3,2] (3,6,2,1,0,10);
array result[2,2];
call mult(mat1, mat2, result);
rc = write_array(&quot;RESULT&quot;, result ) ;
quit;
※ちなみに，関数が成功した場合，戻り値は0
proc fcmp;
array mat1[2,3] (1,2,4,5,8,9);
array mat2[3,2] (3,6,2,1,0,10);
array result[2,2];
call mult(mat1, mat2, result);
rc = write_array(&quot;RESULT&quot;, result ,&quot;COL1&quot;,&quot;COL2&quot;) ;
quit;
※第三引数以降，文字列リストや配列で変数名指定可能


# Page. 6

![Page Image](https://bcdn.docswell.com/page/GE8R5DGLJD.jpg)

関数
説明
READ_ARRAY
データセットをFCMP配列(行列)に変換
第一引数に”入力データセット名”, 第二引数に配列名，第三以降は読み込み変数を指定の場合使用
proc fcmp;
array mat1[2,3] /nosym; ※nosymbols(nosym)は変数や値の割り当
array mat2[3,2] /nosym; て無しで一旦配列定義するための宣言．
array result[2,2];
rc = read_array(&quot;wk1&quot;, mat1) ;
※第二引数以降指定しな場合，データセッ
rc = read_array(&quot;wk2&quot;, mat2 ,&quot;var1&quot;,&quot;var2&quot;) ; トのすべての変数が読み込まれる．部分的
put mat1=;
に配列に入れたい場合は２つめのように変
数名を文字列指定
put mat2=;
call mult(mat1, mat2, result);
rc = write_array(&quot;RESULT&quot;, result ) ;
※⇩データセットの値が配列に格納されている
quit;


# Page. 7

![Page Image](https://bcdn.docswell.com/page/LELKYMGQJR.jpg)

行列操作ルーチン
説明
CALL INV
逆行列を計算
成立要件：入力行列と出力行列は正方で、同じ次元が含まれている必要があります。
※正方(行列)：行要素の数と列要素の数が一致する行列. n×n
proc fcmp;
array mat1[3,3] (3,5,0,9,1,2,8,7,9);
array inv_mat1[3,3] /nosym;
array result[3,3] /nosym;
/*call invでmat1の逆行列をinv_mat1に格納 */
call inv(mat1,inv_mat1);
/*mat1と逆行列の行列積をresultに格納*/
call mult(mat1, inv_mat1, result);
rc = write_array(&quot;MAT1&quot;, mat1 ) ;
rc = write_array(&quot;INV_MAT1&quot;, inv_mat1 ) ;
⇩逆行列
rc = write_array(&quot;RESULT&quot;, result ) ;
quit;
proc print data=MAT1;
run;
proc print data=INV_MAT1;
run;
proc print data=RESULT;
run;
1
-1
↑乗法積が単位行列
-がついてるのは
浮動小数点誤差


# Page. 8

![Page Image](https://bcdn.docswell.com/page/4JM2NYLKJW.jpg)

行列操作ルーチン
説明
CALL TRANSPOSE
行列の転置
proc fcmp;
At などど表記
array mat1[3,2] (3,5,0,9,1,2);
※次元m x nを転置すると次元 n x m になる
array trans_mat1[2,3] /nosym;
/*call transposeでmat1の転置行列をtrans_mat1に格納 */
call transpose(mat1,trans_mat1);
rc = write_array(&quot;MAT1&quot;, mat1 ) ;
rc = write_array(&quot;TRANS_MAT1&quot;, trans_mat1 ) ;
quit;
proc print data=MAT1;
run;
proc print data=TRANS_MAT1;
run;


# Page. 9

![Page Image](https://bcdn.docswell.com/page/PJRMD9K6E9.jpg)

行列操作ルーチン
説明
CALL IDENTITY
入力行列を単位行列に変換
成立要件：入力行列は正方にしてください。
proc fcmp;
array mat1[3,3] (3,5,0,9,1,2,8,7,9);
array mat2[5,5] /nosym;
array mat3[3,3] (1,2,3,4,5,6,7,8,9);
単位行列はその対角
array result[3,3] /nosym;
成分に 1 が並び、他
call identity(mat1);
は全て 0 となる
call identity(mat2);
rc = write_array(&quot;MAT1&quot;, mat1 ) ;
rc = write_array(&quot;MAT2&quot;, mat2 ) ;
rc = write_array(&quot;MAT3&quot;, mat3 ) ;
call mult(mat3, mat1, result);
rc = write_array(&quot;RESULT&quot;, result ) ;
quit;
proc print data=MAT1;
run;
proc print data=MAT2;
=
run;
proc print data=MAT3;
run;
単位行列を書けても値は変わらない．
proc print data=RESULT;
※CALL MULTの頁で説明した計算方法をもう一度みてやってみて
run;
×


# Page. 10

![Page Image](https://bcdn.docswell.com/page/PEXV1QLD7X.jpg)

行列操作ルーチン
説明
CALL DET
行列式を計算
成立要件：入力行列は正方にしてください。
行列Aの行列式は、∣A∣やdet A と表す
proc fcmp;
array mat1[2,2] (2 1 3 4);
call det(mat1, result);
put result=;
quit;
2 1
3 4
2次正方行列の行列式
２*4 -1*3=5
Wikipedia [行列式]より．画像はパブリックドメイン
n次正方行列 A の行列式
detA=
෍ (sgnσ)a1σ(1)a2σ(2)…anσ(n)
σ∈Sn
Sn:置換全体を集めた対称群. sgnσ は置換の符号
↑幾何学的イメージ
𝑎
𝑏
と
が作る平行四辺形の面積
𝑐
𝑑
面積1の正方形をAで線形変換してできる平行四辺形の面積


# Page. 11

![Page Image](https://bcdn.docswell.com/page/3EKQ29LDED.jpg)

行列操作ルーチン
説明
CALL COL
コレスキー分解を計算
成立要件：入力行列と出力行列は正方で、同じ次元が含まれている必要があります。Xは対称正定値で、
Yは下三角行列である必要があります
proc fcmp;
array mat1[3,3] 2 2 3 2 4 2 3 2 6;
array mat2[3,3] / nosym;
行列A を下三角行列Cと，転置したCtの積に分解する操作
call chol(mat1, mat2, 0);
rc = write_array(&quot;MAT1&quot;, mat1 ) ; 三角行列：正方行列が上または下側に値を持ち主対角線
rc = write_array(&quot;MAT2&quot;, mat2 ) ; で反対方向が0になる行列
array trans_mat2[3,3] / nosym;
array validation[3,3] / nosym;
call transpose(mat2, trans_mat2);
call mult(mat2,trans_mat2,validation);
rc = write_array(&quot;TRANS_MAT2&quot;, trans_mat2 ) ;
rc = write_array(&quot;VALIDATIONI&quot;, validation ) ;
run;
proc print data=MAT1;
run;
proc print data=MAT2;
run;
proc print data=TRANS_MAT2;
run;


# Page. 12

![Page Image](https://bcdn.docswell.com/page/L73P4W3PE5.jpg)

必要最低限の機能はそろっているが…
IMLに比べると，操作するための関数(ルーチン)が圧倒的に少ない．
三角行列と簡単に取得したり，ブロック行列を構成したり…
固有値分解や特異値分解，対称行列の対角化など，コードを組むと中々に骨が折れるものが多い…．
固有値分解などはProc Princomp(主成分分析)などのプロシジャ機能を利用したほうが簡易．
固有値と固有ベクトルをpriccompで計算する例
data have(type=cov);
_type_=&#039;COV&#039;;
input m1-m3;
datalines;
0.145373 -0.261593
1.4433602
-0.261593 1.0143407 -0.116446
1.4433602 -0.116446 0.9282068
;
run;
ods select eigenvalues eigenvectors;
proc princomp data=have(TYPE=COV) cov;
run;
ただし，多変量解析や機械学習などが，複雑な行列操作を伴うことが多く，治験分野の解析やシミュ
レーションであれば…，FCMPでも やってやれんことはない？？？


# Page. 13

![Page Image](https://bcdn.docswell.com/page/87DLQK43EG.jpg)

練習１ 連立方程式の解
中学生の数学問題から．3連連立方程式をFCMPで解いてみる．実務ではないけど行列で値を求める練習
𝑥 − 2𝑦 + 8𝑧 = 7
ቐ2𝑥 + 3𝑦 − 5𝑧 = 14
3x + 4y +19z = −29
1 −2 8
2 3 −5
3 4 19
𝑥
7
𝑦 = 14
𝑧
−29
𝑥
1 −2 8
𝑦 = 2 3 −5
𝑧
3 4 19
-1
7
14
−29
proc fcmp;
array mat1 [3, 3](1 -2 8 2 3 -5 3 4 19);
array r
[3, 1] (7 14 -29);
array _inv [3, 3] / nosym;
array result [3,1]
/ nosym;
call inv(mat1, _inv); ;
call mult(_inv, r, result); ;
rc = write_array(&quot;RESULT&quot;, result) ;
run;
proc print data=RESULT;
run;
A. x = 11, y = -6, z = -2


# Page. 14

![Page Image](https://bcdn.docswell.com/page/VJP2LKMPJ8.jpg)

練習2 重回帰分析
SASHELP.CLASSデータセットを使い．体重を応答変数，身長と年齢を説明変数とする重回帰モデルで
各説明変数の偏回帰係数を推定する．
偏回帰係数ベクトルをβは以下の式で求まる
𝑏0
β= … =ሺ𝑋 𝑡 𝑋)−1 𝑋 𝑡 𝑌
𝑏𝑝−1
proc sql noprint;
/*obs数をマクロ変数に格納*/
select count(*) into:obs from
sashelp.class;
quit;
data x;
set sashelp.class ;
x1=1; /*切片用変数、全obs共通で1を格納*/
x2=HEIGHT; /*身長*/
x3=AGE; /*年齢*/
keep x: ;
run;
data y;
set sashelp.class;
y = WEIGHT; /*体重*/
keep y;
run;


# Page. 15

![Page Image](https://bcdn.docswell.com/page/2EVDQV9VEQ.jpg)

練習2 重回帰分析
proc fcmp;
array x[&amp;obs.,3]/nosym;
𝑏0
array y[&amp;obs.,1]/nosym;
𝑡
−1 𝑡
…
ሺ
β=
=
𝑋
𝑋)
𝑋 𝑌
array beta[3,1]/nosym;
𝑏𝑝−1
/*説明変数と応答変数それぞれ読み込み*/
rc=read_array(&#039;x&#039;,x);
rc=read_array(&#039;y&#039;,y);
proc reg data=sashelp.class;
proc print data=beta;
/*①xの転置*/
model WEIGHT = HEIGHT AGE;
run;
array tran_x[3,&amp;obs.]/nosym;
quit;
call transpose (x,tran_x);
/*②xと転置xの積*/
array x_tran_x[3,3]/nosym;
call mult(tran_x,x,x_tran_x);
/*③(xと転置xの積)の逆行列*/
array x_tran_x_inv[3,3]/nosym;
call inv(x_tran_x,x_tran_x_inv);
/*④3と転置xの積*/
array x_tran_x_inv_tran_x[3,&amp;obs.]/nosym;
call mult(x_tran_x_inv,tran_x,x_tran_x_inv_tran_x);
/*⑤4とyの積をβに格納*/
call mult(x_tran_x_inv_tran_x, y, beta);
rc=write_array(&#039;beta&#039;,beta);
run;


# Page. 16

![Page Image](https://bcdn.docswell.com/page/57GYWLZ17L.jpg)

練習3 対称行列の対角化
data m;
input c1-c5;
cards;
10 11 12 13 14
11 20 11 11 11
12 11 30 15 16
13 11 15 40 25
14 11 16 25 50
;
run;
data cov(type = cov);
length _type_ _name_ $20;
array ar[*] c1-c5;
set m end = eof;
_type_ = &#039;cov&#039;;
_name_ = vname(ar[_n_]);
output;
if eof then do;
_type_ = &#039;N&#039;;
_name_ = &#039;&#039;;
call stdize(&#039;mult=&#039;, 0, &#039;add=&#039;, 100, of
ar[*]);
output;
end;
run;
ods output eigenvalues = eigenvalues eigenvectors =
eigenvectors;
proc princomp data = cov cov;
run;
ods output close;
proc fcmp;
array _V[5, 5] / nosym;
array _U[5, 5] / nosym;
array _Up[5, 5] / nosym;
array _VU[5, 5] / nosym;
array _UpVU[5, 5] / nosym;
rc = read_array(&#039;m&#039;, _V, &#039;c1&#039;, &#039;c2&#039;, &#039;c3&#039;, &#039;c4&#039;, &#039;c5&#039;);
rc = read_array(&#039;eigenvectors&#039;, _U, &#039;prin1&#039;, &#039;prin2&#039;, &#039;prin3&#039;,
&#039;prin4&#039;, &#039;prin5&#039;);
call transpose(_U, _Up);
call mult(_V, _U, _VU);
call mult(_Up, _VU, _UpVU);
rc = write_array(&#039;out&#039;, _UpVU, &#039;_1&#039;, &#039;_2&#039;, &#039;_3&#039;, &#039;_4&#039;, &#039;_5&#039;);
run;


# Page. 17

![Page Image](https://bcdn.docswell.com/page/4EQM3YLNJP.jpg)

所感
proc fcmp内で提供されている行列操作ルーチンを利用することで，基本的な行列計算は十分
可能である．
IMLが利用できない環境で，簡易な行列計算が必要な場合，有用な手段と考える．
また，今回，いいサンプルを提示することができなかったが，FCMPはもともとユーザ定義関
数を定義可能なプロシジャのため，行列計算を含むユーザー定義関数の定義などはIMLにはな
い長所となりえる．
IMLに不足している機能を，ユーザー定義関数で補うことも可能であろう．
一方で，可変的な引数と戻り値をもったユーザー定義関数やユーザー定義サブルーチンを作成
する難易度はかなり高く，悩ましいところである．
昨今では，シミュレーションや，複雑な行列計算を含む処理にRを利用することも増えている
と思われる．
しかし，常に選択肢は多いに越したことがない．
行列処理の一つの手段としてFCMPプロシジャの利用は十分一考に値すると考えている


# Page. 18

![Page Image](https://bcdn.docswell.com/page/KJ4V1WD3E1.jpg)

おまけ
Proc Luaで行列操作関数を一通り


# Page. 19

![Page Image](https://bcdn.docswell.com/page/LE13GYZZJG.jpg)

LuaはPythonとかと似た感じ（暴言）
やろうと思えば行列操作の関数群，パッケージを作ること
は比較的簡単
⇒
SASの何倍も情報があるので，LLMでの生成精度がよくて
それをproc luaの中に貼ればおおむね動いてしまう
以下はあくまで参考で．間違ってたらごめんなさい


# Page. 20

![Page Image](https://bcdn.docswell.com/page/GEWPKG96J2.jpg)

proc lua;
submit;
-- 行列操作ライブラリ
Matrix = {}
-- 新しい行列を作成
function Matrix.new(rows, cols, value)
local mat = {}
value = value or 0
for i = 1, rows do
mat[i] = {}
for j = 1, cols do
mat[i][j] = value
end
end
return mat
end
-- 行列の加算
function Matrix.add(a, b)
local rows = #a
local cols = #a[1]
if rows ~= #b or cols ~= #b[1] then
error(&quot;Matrix dimensions must match for addition&quot;
)
end
local result = Matrix.new(rows, cols)
for i = 1, rows do
for j = 1, cols do
result[i][j] = a[i][j] + b[i][j]
end
end
return result
end


# Page. 21

![Page Image](https://bcdn.docswell.com/page/47ZKZL9R73.jpg)

-- 行列の乗算
function Matrix.multiply(a, b)
-- 行列の減算
local rows_a = #a
function Matrix.subtract(a, b)
local cols_a = #a[1]
local rows = #a
local rows_b = #b
local cols = #a[1]
local cols_b = #b[1]
if rows ~= #b or cols ~= #b[1] then
error(&quot;Matrix dimensions must match for subtracti
if cols_a ~= rows_b then
on&quot;)
error(&quot;Matrix dimensions incompatible for multiplic
end
ation&quot;)
end
local result = Matrix.new(rows, cols)
for i = 1, rows do
local result = Matrix.new(rows_a, cols_b)
for j = 1, cols do
for i = 1, rows_a do
result[i][j] = a[i][j] - b[i][j]
for j = 1, cols_b do
end
local sum = 0
end
for k = 1, cols_a do
return result
sum = sum + a[i][k] * b[k][j]
end
end
result[i][j] = sum
end
end
return result
end


# Page. 22

![Page Image](https://bcdn.docswell.com/page/YJ6QZWK1JV.jpg)

-- スカラー倍
function Matrix.scalarMultiply(mat, scalar)
local rows = #mat
local cols = #mat[1]
local result = Matrix.new(rows, cols)
-- 行列の転置
function Matrix.transpose(mat)
local rows = #mat
local cols = #mat[1]
local result = Matrix.new(cols, rows)
for i = 1, rows do
for j = 1, cols do
result[i][j] = mat[i][j] * scalar
end
end
return result
end
for i = 1, rows do
for j = 1, cols do
result[j][i] = mat[i][j]
end
end
return result
end


# Page. 23

![Page Image](https://bcdn.docswell.com/page/GJ5RWMPL74.jpg)

-- アダマール積（要素ごとの積）
function Matrix.hadamard(a, b)
local rows = #a
local cols = #a[1]
if rows ~= #b or cols ~= #b[1] then
error(&quot;Matrix dimensions must match for Hadamard produ
ct&quot;)
end
local result = Matrix.new(rows, cols)
for i = 1, rows do
for j = 1, cols do
result[i][j] = a[i][j] * b[i][j]
end
end
return result
end
-- 単位行列の生成
function Matrix.identity(n)
local result = Matrix.new(n, n)
for i = 1, n do
result[i][i] = 1
end
return result
end
-- 逆行列の計算（ガウス・ジョルダン法を使用）
function Matrix.inverse(mat)
local n = #mat
if n ~= #mat[1] then
error(&quot;Matrix must be square for inverse&quot;)
end
-- 拡張行列を作成 [A|I]
local aug = Matrix.new(n, 2*n)
for i = 1, n do
for j = 1, n do
aug[i][j] = mat[i][j]
end
aug[i][i + n] = 1 -- 単位行列部分
end
-- ガウス・ジョルダン法
for i = 1, n do
-- ピボットの選択（簡易版）
local pivot = aug[i][i]
if pivot == 0 then
error(&quot;Matrix is singular or nearly singular&quot;)
end
-- ピボット行を1にする
for j = 1, 2*n do
aug[i][j] = aug[i][j] / pivot
end
-- 他の行を0にする
for k = 1, n do
if k ~= i then
local factor = aug[k][i]
for j = 1, 2*n do
aug[k][j] = aug[k][j] - factor * aug[i][j]
end
end
end
end


# Page. 24

![Page Image](https://bcdn.docswell.com/page/9E2DQ963ER.jpg)

-- 対角行列の生成
function Matrix.diagonal(values)
local n = #values
local result = Matrix.new(n, n)
for i = 1, n do
result[i][i] = values[i]
end
return result
end
-- 使用例
local a = Matrix.new(2, 2)
a[1][1] = 1; a[1][2] = 2
a[2][1] = 3; a[2][2] = 4
local b = Matrix.new(2, 2)
b[1][1] = 5; b[1][2] = 6
b[2][1] = 7; b[2][2] = 8
print(&quot;Matrix A:&quot;)
Matrix.print(a)
print(&quot;\nMatrix B:&quot;)
Matrix.print(b)
print(&quot;\nA + B:&quot;)
local sum = Matrix.add(a, b)
Matrix.print(sum)
print(&quot;\nA * B:&quot;)
local product = Matrix.multiply(a, b)
Matrix.print(product)
print(&quot;\n2 * A:&quot;)
local scaled = Matrix.scalarMultiply(a, 2)
Matrix.print(scaled)
print(&quot;\nTranspose of A:&quot;)
local transposed = Matrix.transpose(a)
Matrix.print(transposed)
- アダマール積
print(&quot;\nHadamard Product:&quot;)
local had = Matrix.hadamard(a, b)
Matrix.print(had)
-- 単位行列
print(&quot;\nIdentity Matrix (3x3):&quot;)
local I = Matrix.identity(3)
Matrix.print(I)
-- 逆行列
print(&quot;\nOriginal Matrix A:&quot;)
Matrix.print(a)
print(&quot;\nInverse of A:&quot;)
local inv = Matrix.inverse(a)
Matrix.print(inv)
print(&quot;\nA * A^-1 (should be identity):&quot;)
local check = Matrix.multiply(a, inv)
Matrix.print(check)
endsubmit;
quit;


# Page. 25

![Page Image](https://bcdn.docswell.com/page/D7Y2W4987M.jpg)



# Page. 26

![Page Image](https://bcdn.docswell.com/page/VEN89YL9J8.jpg)

-- SASデータセットをLua行列に変換する関数
function Matrix.sasDatasetToMatrix(datasetName)
-- データセットを開く
local dsid = sas.open(datasetName)
if not dsid then
error(&quot;Cannot open dataset: &quot; .. datasetName)
end
-- 変数のリストを取得
local varNames = {}
for var in sas.vars(dsid) do
table.insert(varNames, var.name)
end
print(table.tostring(varNames))
-- 行列を格納するテーブル
local matrix = {}
-- データセットの各行を読み込む
local rowNum = 1
while sas.next(dsid) do
matrix[rowNum] = {}
for i, varName in ipairs(varNames) do
matrix[rowNum][i] = sas.get_value(dsid, varName)
end
rowNum = rowNum + 1
end
-- データセットを閉じる
sas.close(dsid)
return matrix, varNames
end
逆に，行列をSASデータセットにする場合は，
Proc Lua組み込み関数のsas.write_dsで一応事足りるので
そっちつかってください


# Page. 27

![Page Image](https://bcdn.docswell.com/page/Y79K2P4373.jpg)

さらにおまけ
ガウス・ニュートン法で方程式
（単変数関数）の解を求めることを
FCMP.SOLVE関数で実装する方法


# Page. 28

![Page Image](https://bcdn.docswell.com/page/G78R5DQLED.jpg)

1
ｙ=
𝑥
y=20になるxを求める例
data wk1;
do x =0 to 0.01 by 0.0001;
y=1/sqrt(x);
output;
end;
run;
proc sgplot data=wk1;
refline 20;
series x=x y=y;
run;


# Page. 29

![Page Image](https://bcdn.docswell.com/page/L7LKYMXQER.jpg)

proc fcmp;
/*解きたい式を関数定義*/
function f1(x);
y = 1/sqrt(x);
return(y);
endsub;
1
ｙ=
𝑥
/*solve関数 */
x=solve(&quot;f1&quot;,{.}, 20 , .);
put x= ;
run;
sovle関数の第一引数に与えます．第二引数の{.} はガウス・ニュートン法で 解を求める際の初期
値や，反復数，どれぐらいまでの近似精度にするかなどのパラメータ設定を配列で渡す部分で，空
配列を指定するとすべてデフォルト値を採用します


# Page. 30

![Page Image](https://bcdn.docswell.com/page/4EM2NY3KEW.jpg)

proc fcmp outlib=work.common.solve ;
/*解きたい式を定義*/
function f1(x);
y = 1/sqrt(x);
return(y);
endsub;
/*yを指定して、xを求めて戻すルーチン*/
subroutine solver(y,x);
outargs x ;
x=solve(&quot;f1&quot;,{.}, y , .);
endsub ;
run;
options cmplib=work.common;
data wk1;
call missing(x);
call solver(20,x);
run;
proc print data=wk1;
run;


# Page. 31

![Page Image](https://bcdn.docswell.com/page/PERMD9Q679.jpg)

proc fcmp;
/*解きたい式を定義*/
function f1(x);
y = 1/sqrt(x);
return(y);
endsub;
y=20;
/*設定パラメータ */
array opts[5] initial_value
absolute_criterion
relative_criterion
maximum_iterations
solve_status;
initial_value=0.001;
absolute_criterion=1.0e-12;
relative_criterion=1.0e-6;
maximum_iterations=100;
solve_status=.;
initial_valueは初期値で，これがあまりにも離れた値だ
と収束しない可能性があります
また複数解をもつ場合に，どの解が選択されるかは初
期値の影響をうけます
absolute_criterionは目標としてる値（今回20）との誤
差をどの程度で，一致とみなすか
relative_criterionは計算された値の変化率が解であろう
付近極小である場合，どこまでの安定で収束とみなす
かです（多分）
maximum_iterationsは最大反復数
solve_statusは，欠損で指定しておけば あとで結果
に応じて 0なら成功，それ以外なら
x=solve(&quot;f1&quot;
,opts
, y
, .);
put x= solve_status=;
run;
quit;
この引数は対象の関数に対して，引数を与えるときに使う


# Page. 32

![Page Image](https://bcdn.docswell.com/page/P7XV1Q5DJX.jpg)

proc fcmp outlib=work.common.solve ;
/*解きたい式を定義*/
function f1(x);
return(x**3 -6*x**2 + 11*x -6);
endsub;
/*yを指定して、xを求めて戻すルーチン*/
subroutine
solver(y,x,initial_value,absolute_criterion,relative_cr
iterion,maximum_iterations,solve_status);
/*設定パラメータ */
array opts[5] initial_value absolute_criterion
relative_criterion maximum_iterations solve_status;
outargs x ,solve_status,y ;
x=solve(&quot;f1&quot;
,opts
, 0
, .);
put x= solve_status=;
endsub ;
run;
options cmplib=work.common;


# Page. 33

![Page Image](https://bcdn.docswell.com/page/37KQ296D7D.jpg)

options cmplib=work.common;
data wk1;
call missing(of x);
absolute_criterion=1.0e-12;
relative_criterion=1.0e-6;
maximum_iterations=100;
solve_status=.;
y=0;
グラフから，おそらくｘ＝1,2,3あたりに解がありそうなので
初期値(initial_value)を，0.001，1.9, 2.9として回してみると
initial_value=0.001;
call
solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterat
ions,solve_status);
output;
initial_value=1.9;
call
solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterat
ions,solve_status);
output;
initial_value=2.9;
call
solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterat
ions,solve_status);
output;
run;
proc print data=wk1;
run;


# Page. 34

![Page Image](https://bcdn.docswell.com/page/LJ3P4W6P75.jpg)

やはり ｘ=1,2,3が解でした．solve_status=0は収束してい
ることを示します
ガウス・ニュートン法のアルゴリズム自体の説明は省きます
関数の最大・最小値を見出すニュートン法の改良で，非線形最小二乗法を解く手法の一つです．
全ての状況で解を求めるのに使えるわけではないですが，適用可能な場合もあるので
選択肢の一つとしていいですね


