行列の関数とクラス

2017-08-05

JavaScript における行列

前回は JavaScript で多項式のクラスを考えた。 今度は行列のクラスを考えようとした。しかし、行列はある程度 JavaScript の 2 次元配列で操作できる。 まずは関数を実装することを先にして、クラスは後回しにする。

ここで、行列とは、`A = [(a_(11), a_(12), a_(13) ) , ( a_(21), a_(22), a_(23) )] ` のようなものをいう。このような行列を操作する関数を作る。

行列を二次元配列（多重配列）で表しておけば、JavaScript の配列のメソッドが使える。

演算

以下の説明では、JavaScript による mathjs というライブラリで用意されている関数を使った結果を合わせて説明する。

ベクトルの内積

i と j を固定してa[i] と a[j] をそれぞれベクトルとみなしたときの内積を作る。次は i = 1, j = 2 とした場合だ。 

const a=[[1,2,3],[4,5,6],[7,8,9]]
const f = a[1].reduce((is, as, i) => is + as * a[2][i], 0)
console.log(f) // => 122 (=28+40+54)

内積を作る関数に ip (エス : skalara produto, punkta produto, interna produto) という名前をつけよう。 

const a=[[1,2,3],[4,5,6],[7,8,9]]
const ip = (v, w) => v.reduce((is, as, i) => is + as * w[i], 0);
console.log(ip(a[1], a[2])) // => 122 (=28+40+54)

次に、内積を発展させよう。2 次元配列 a に対し、1 次元配列 a[0] を固定し、 a[0] とそれぞれ配列 a[0], a[1], a[2] の内積をそれぞれ求める。 

const a=[[1,2,3,4],[5,6,7,8],[9,10,11,12]]
const g = a.map(v => ip(v, a[0]));
console.log(g) // => [30, 70, 110] ([1*1+2*2+3*3+4*4,1*5+2*6+3*7+4*8,1*9+2*10+3*11+4*12])

mathjs で内積を計算することもできる。

加減算

`C = A + B , D = A - B` を計算する。 `A = [[1,4,2], [8,5,7]]`,`B = [[8,5,7], [1,4,2]]` としたとき、 `C = [[9,9,9], [9,9,9]], D = [[-7,-1,-5],[7,1,5]]` である。

mathjs では次のようになる。

math.add([[1,4,2], [8,5,7]],[[8,5,7], [1,4,2]]) // returns [[9,9,9],[9,9,9]]

math.subtract([1,4,2], [8,5,7]) // returns [[9,9,9],[9,9,9]]

乗算

`G = EF` を計算する。 `E = [[1,-4,0,0], [1,0,5,-3],[0,1,7,2],[9,-3,0,0]]`, `F = [[4,0,0], [0,4,0],[0,0,-3],[0,0,-1]]` としたとき、 `G = [[4,-16,0], [4,0,-12],[0,4,-23],[36,-12,0]]` である。

mathjs では次のようになる。

math.multiply([[1,-4,0,0], [1,0,5,-3],[0,1,7,2],[9,-3,0,0]], [[4,0,0], [0,4,0],[0,0,-3],[0,0,-1]]) //
[[4,-16,0], [4,0,-12],[0,4,-23],[36,-12,0]]

小行列

`i` 行 `j` 列を除いた行列を小行列という。行列 `A` の `i` 行 `j` 列を除いた行列を `A(i, j)` と書く。 上記の `E` に対して、`H= E(0,0)` とおくと `H` は次のようになる。
`H = [[0,5,-3],[1,7,2],[-3,0,0]]`。

mathjs による小行列を求めるのは少し厄介だ。削除する行や列を指定する方法は後回しにして、 まず残す行や列のインデックスを指定して小行列を求める関数を説明する。

math.subset([[1,2,3],[4,5,6],[7,8,9]], math.index([0,1,3], [0,2])) // returns [[1,3],[4,6],[10,12]]

これを利用して改めて削除する行や列の番号を指定する方法を考える。回りくどいが次のように subset と index, concat, range を組み合わせればできる。 スマートなやり方は私にはわからない。 

const a = [[1,2,3],[4,5,6],[7,8,9],[10,11,12]]
	const i = 2 // 3 行目を消す
	const j = 1 // 2 列目を消す
	const r_index = math.concat(math.range(0, i),math.range(i + 1, 4)) // 全部で 4 行
	const c_index = math.concat(math.range(0, j),math.range(j + 1, 3)) // 全部で 3 列
	math.subset(a, math.index(r_index, c_index)) // [[1,3],[4,6],[10,12]]

行列式

行列式の定義にはいろいろあるが、省略する。

mathjs では次のようになる。

math.det([[2, -1, 3, -2], [1, 7, 1, -1], [3, 5, -5, 3], [4, -3, 2, -1]]) // 42

以下は、数学で使われるベクトルや行列の演算を例として、上記で紹介した各種 JavaScript の関数やメソッドの使い方を確認する。

行方向・列方向への加算

行列では、行方向への加算や列方向への加算を考えることができる。まず、行方向への加算について考える。

b[i]=a[i][0]+a[i][1]+a[i][2] となる配列 b[i] (i=0, 1, 2)を作る

const a=[[1,2,3],[4,5,6],[7,8,9]]
const b = a.map(h => h.reduce((is, as) => is + as, 0) ) // is = prev, as = current 
console.log(b) // => [6, 15, 24]

次に列方向への加算を考える。

d[i]=a[0][i]+a[1][i]+a[2][i] となる配列 d[i] (i=0, 1, 2)を作る

const a=[[1,2,3],[4,5,6],[7,8,9]]
const d = a.reduce((is, as) => is.map((v,i) => v+as[i]), a[0].slice().fill(0) );
console.log(d) // => [12, 15, 18]

この d を作る例は私には難しい。 [1,2,3]+[4,5,6] ができるようなコールバック関数を reduce に与えなければならないので、 このように map を使っている。reduce で使える引数を is, as で与えているが、 これはエスペラントの過去時制を表わす接尾辞 -is と現在時制を表わす接尾辞 -as である。 英語ならばそれぞれ previousValue, currentValue を使うだろう。 もう一つ注意すべき点は、reduce に与える初期値を a[0].slice().fill(0) としたことである。 最初 a[0].fill(0) としていたが、これではもとの a[0]の値を上書きしてしまうので、 誤った値 [11, 13, 15] になってしまう。そこで slice を使ってコピーを作るようにした。 2022-04-06

行方向や列方向の加算は、画像処理ではプロジェクションとも呼ばれる。

列の抜き出し

c[i]=a[0][i] となる配列 c[i] (i=0, 1, 2)を作る

const a=[[1,2,3],[4,5,6],[7,8,9]]
const c = a.map(v => v[0])   
console.log(c) // => [1, 4, 7]

上記は 1 列目を取り出す方法だ。同様に 2 列目、 3 列目もできる。これがわかると、行列の転置の方法もわかるはずだ。

なお、プログラムを使わずに、mathjs で提供されているライブラリを使う列の抜き出しもある。

行列の転置

a[i][j] の転置行列を作る。 すなわち e[j][j]=a[i][j] となる配列 e[i] を作る。ここでは、i, j = 0, 1, 2 とする。なお、以下の方法では、行数と列数が異なる場合でも問題ない。

const a=[[1,2,3],[4,5,6],[7,8,9]]
const e = a[0].map((_, j) => a.map(h => h[j]))
console.log(e) // => [[1,4,7],[2,5,8],[3,6,9]]

この例も私には難しい。最初の map メソッドの引数 (_, j) は、第 1 の引数は捨て、 第2の引数 j のみを使うことを意味する。j は列が j 番目であることを表している。 第 j 列に対して再度 a.map で横（行）horizontalo ごとに h[j] を呼び出してベクトルを作る。 このベクトルを作ることについては、列の抜き出しも参照されたい。 なお、行列の転置を関数とすることもできる。仮に trans という名前をつけよう。 

const a=[[1,2,3],[4,5,6],[7,8,9]]
const trans = m => m[0].map((_, v) => m.map(h => h[v]));
console.log(trans(e)) // => [[1,4,7],[2,5,8],[3,6,9]]

三角行列の転置

a が三角行列である場合の転置行列を作る。まず a が下三角行列である場合を考える。つまり、下三角から上三角への変換である。

const a=[[1],[2,3],[4,5,6]]
const e = a.map((_, i) => a.slice(i).map(h => h[i]))
console.log(e) // => [[1,2,4],[3,5],[6]

次に a が上三角行列である場合を考える。

const a=[[1,2,3],[4,5],[6]]
const e = a.map((_, i) => a.slice(0, i + 1).map((h, j) => h[i - j]))
console.log(e) // => [[1],[2,4],[3,5,6]

この上三角から下三角への変換は、その逆に比べて難しい。これらの変換は、a.slice(i) や a.slice(0, i + 1) が肝である。 前者は i が増えるにつれて a を上から削る量を多くしている。一方後者は i が増えるにつれて a を下から削る量を少なくしている。 また、前者では不要だった 内側の map の j が、後者では必要となっている。これについて細かな説明は省略する。

行列の積

さらに発展させて、行列の積を作ることができる。

const a=[[1,2,3,4],[5,6,7,8],[9,10,11,12]]
const b=[[3,4],[5,6],[7,8],[9,10]];
const h = a.map(v => trans(b).map(w => ip(v, w)));
console.log(h) // => [[70, 80],[166,192],[262,304]]

特に、行列が 1 つの場合、a を行列として a と trans(a) の積を作る。 

const a=[[1,2,3,4],[5,6,7,8],[9,10,11,12]]
const k = a.map(v => a.map(w => ip(v, w)));
console.log(k) // => [[30, 70, 110],[70, 174, 278],[110, 270, 446]]

この結果を見るとわかる通り、計算には無駄がある。積は対称行列だから、計算は上三角部分だけでいい。 a[i]とa[j]の内積を a[i]*a[j] のように書くとき、次のような二次元配列が作れるか。
[[a[0]*a[0],a[0]*a[1],a[0]*a[2]],[a[1]*a[1],a[1]*a[2]],[a[2]*a[2]]]
妥協して次の書き方にした。 

const a=[[1,2,3,4],[5,6,7,8],[9,10,11,12]]
const l = a.map(v, i) => a.map((w, j) => (i <= j ? ip(v, w) : undefined));
console.log(l) // => [[30, 70, 110],[, 174, 278],[,, 446]]

このプログラムを書いたとき「もう少しスマートな書き方があるはずだ」と記した後１年以上考えた。 その結果、上三角行列を得るには次のように書けばいいことがわかった。ここで ip(w, u) は配列 w と配列 u の内積をとる関数である。 (2025-01-08)


const a=[[1,2,3,4],[5,6,7,8],[9,10,11,12]]
const l = a.map((u, i, v) => v.slice(i).map(w => ip(w, u)))
console.log(l) // => [[30, 70, 110],[174, 278],[446]]

下三角行列を得るには次のように書けばいい。(2025-01-08)


const a=[[1,2,3,4],[5,6,7,8],[9,10,11,12]]
const l = a.map((u, i, v) => v.slice(0, i + 1).map(w => ip(w, u)))
console.log(l) // => [[30],[70, 174],[110, 278, 446]]

拡大行列

二つの与えられた行列の列を連結させて得られる行列を拡大行列という。 拡大行列は、二つの行列に対し同じ行基本変形を施すことを目的として構成される。 下記の例では列の数を合わせるために、b = [[4],[3],[1]] としている。 b = [4,3,1] ではないことに注意。

const a=[[1,3,2],[2,0,1],[5,2,2]], b = [[4],[3],[1]];
const c = a.map((v, i) => v.concat(b[i]));
console.log(c); // [[1,3,2,4],[2,0,1,3],[5,2,2,1]]

ただ、この場合は、b = [4,3,1] でも同じ結果が得られる。

const a=[[1,3,2],[2,0,1],[5,2,2]], b = [4,3,1];
const c = a.map((v, i) => v.concat(b[i]));
console.log(c); // [[1,3,2,4],[2,0,1,3],[5,2,2,1]]

行方向に連結させるときは注意が必要である。次はうまくいかない。

const a=[[1,3,2],[2,0,1],[5,2,2]], b = [4,3,1];
const c = a.concat(b);
console.log(c); // [[1, 3, 2], [2, 0, 1], [5, 2, 2], 4, 3, 1]

次もうまくいかない。

const a=[[1,3,2],[2,0,1],[5,2,2]], b = [[4],[3],[1]];
const c = a.concat(b);
console.log(c); // [[1, 3, 2], [2, 0, 1], [5, 2, 2], 4, 3, 1]

次はうまくいく。

const a=[[1,3,2],[2,0,1],[5,2,2]], b = [[4,3,1]];
const c = a.concat(b);
console.log(c); // [[1, 3, 2], [2, 0, 1], [5, 2, 2], [4,3,1]]

次のようにする手もある。

const a=[[1,3,2],[2,0,1],[5,2,2]], b = [4,3,1];
const c = a.concat(Array.of(b));
console.log(c); // [[1, 3, 2], [2, 0, 1], [5, 2, 2], [4,3,1]]

前進消去と後退代入

連立一次方程式を解くときに、行列を LU 分解したりコレスキー分解したあとで解を求める方法が前進消去と後退代入である。 奥村晴彦「C 言語による標準アルゴリズム事典」のp.384 「LU 分解」の例を挙げる。次の例は `A = LU` となっている。

`A = [[2,5,7],[4,13,20],[8,29,50]], L = [[1,0,0],[2,1,0],[4,3,1]], U = [[2,5,7],[0,3,6],[0,0,4]]`

`Ax = b` の形の連立一次方程式を解くには、まず `(LU)x = b` から `L(Ux) = b` と変形する。`Ux = y` とおくと `Ly = b` であり、`L` の形から `y` を求めることができる。その後で `Ux = y` から `x` を求めることができる。

前進消去

まず、`Ly = b` から `y` を求めてみよう。これには前進消去という方法を使う。

`L[0][0]*y[0] = b[0]`
`L[1][0]*y[0]+L[1][1]*y[1] = b[1]`,
`L[2][0]*y[0]+L[2][1]*y[1]+L[2][2]*y[2]=b[2]`
であるから、`L[i][i]=1`と合わせて、次の式が得られる。

`y[0] = b[0], y[1] = b[1] - L[1][0]*y[0], y[2] = b[2] - L[2][0]y[0] - L[2][1]y[1]`

プログラムは次のようになるだろう。


const n = b.length;
for (let i = 0; i < n; i++) {
  y[i] = b[i];
  for (let j = 0; j < i; j++) {
    y[i] -= L[i][j] * y[j];
  }
}

`b = [23,58,132]` として計算してみると、`y=[23,12,4]` が得られる。

この方法は確実だが、インデックスを使ってループを回している。他の方法はどうだろうか。j のループは reduce が使えそうだ。 i のループはひとまず手で展開すれば、次のようなプログラムが考えられる。 

			y = b.slice();
			y[0] = L[0].slice(0, 0-n).reduce((acc, val, j) => acc - val * y[j], b[0])
			y[1] = L[1].slice(0, 1-n).reduce((acc, val, j) => acc - val * y[j], b[1])
			y[2] = L[2].slice(0, 2-n).reduce((acc, val, j) => acc - val * y[j], b[2])

この方法でも `y=[23,12,4]` が得られる。それでは、i のループのかわりに map を使ってみよう。

  let y = b.slice();
  y = L.map((l, i) => l.slice(0, i - n).reduce((acc, val, j) => acc - val * y[j], b[i]) );

この結果は `y=[23, 12, -134]` となり、誤った値となっている。`y` の扱い方に問題があるのだろう。 map のかわりに forEach を使ってみた。次はどうだろうか。

	let y = b.slice();
	L.forEach((l, i) => { y[i] = l.slice(0, i - n).reduce((acc, val, j) => acc - val * y[j], b[i]) })

for ループを使うのとほとんど変わりないが、`y=[23,12,4]` という正しい値が得られている。

後退代入

なお、上三角行列 `U` に対して `Ux = y` から `x` を求める方法は後退代入と呼ばれる。こちらについても見てみよう。

`U[0][0]*x[0] + U[0][1]*x[1] + U[0][2] * x[2] = y[0]`,
` U[1][1]*x[1] + U[1][2] * x[2] = y[1]`,
` U[2][2] * x[2] = y[2]`
であるから、次の式が得られる。

`x[2] = (y[2]) / (U[2][2]), x[1] = (y[1] - U[1][2]*x[2]) / (U[1][1]), x[0] = (y[0] - U[0][2]*x[2] - U[0][1]*x[1]) / (U[0][0])`

プログラムは次のようになるだろう。


const n = b.length - 1;
for (let i = n; i ≥ 0; i--) {
  x[i] = y[i];
  for (let j = n; j > i ; j--) {
    x[i] -= U[i][j] * x[j];
  }
  x[i] /= U[i][i];
}

後退代入と reduceRight メソッド

前進消去で reduce メソッドを使う例を紹介した。そこで、後退代入で recudeRight メソッドを使う例を紹介する。

reduceRightの適用例として、ガウスの消去法や LU 分解、QR 分解などで頻出する後退代入を取り上げよう。具体的には `[[1,2,-1],[0,1,3],[0,0,1]] [[x_1],[x_2],[x_3]] = [[4],[5],[1]]` のような連立一次方程式を解く方法である。解は、
`x_3 = 1 / 1 = 1`
`x_2 = (5 - 3*x_3) / 1 = 2`
`x_1 = (4 - (2*x_2 - x_3)) / 1 = 1`
これを JavaScript で組むと次のようになる。`U` は `n` 次上三角行列を表す配列 （第1行が長さ `n` の配列、第 2 行が長さ `n-1` の配列、`cdots`、第 `n` 行が長さ `1` の配列）、 `x` は未知の `n` 次配列、`b` は既知の `n` 次配列である。上記の場合を例にとると次のようになる。 

const ip = (v, w) => a.reduce((akum, val, i) => akum + val * w[i], 0) // v と w の内積
const U = [[1,2,-1],[1,3],[1]];
const y = [4,5,1];
const x = U.reduceRight((akum, val, i) => {
	return [...akum, (y[i] - ip(akum, val.reverse())) / val[0] ]
}, []).reverse()

ほとんど判じ物である。以下注釈を付ける。なお、内積をとる関数 ip については、内積の項を参照。

二次元配列 `U` は上三角行列を JavaScript で表現したものである。 後退代入なので、配列の最後の要素から取り扱うのが素直だろう。そこで reduceRight を使う。 `x` の最終要素から決めていくわけだが、これは `u_(33) x_3 = y_3` から `x_3 = y_3 // u_(33)` である。 なお、`y_3` は y[i] で表現されている。なぜなら、 a.reduce のcallbackFunction の引数群 akum, val, i のうち i は 配列を array とすると array.length - 1 から始まるからである（初期値が指定されているため）。 それを val[0] で割っている。val は reduce で得られる最初の要素であり、この最初の要素は長さ 1 の配列だから、 対角要素である。よってこれが最初の解の最終要素 `x_3` となり、これが akum に蓄えられる。この結果、reduceRight のループが１回まわり、akum の内容が [] から [`x_3`] となる。

そして２回めのループに入る。return が返す配列の意味であるが、 最初の ...akum は、配列 akum を新たな配列に挿入するために展開していることを表す。そして次の要素は、y[i] から、akum （ここでは長さ1の配列 [`x_3`]） と val.reverse （ここでは長さ 2 の配列[3, 1]）の内積を引いた数を計算し、さらにこの引いた数を対角要素 val[0] で割った値となる。 ここで、内積に関して注意する。 第１引数と第２引数の長さが異なるが、ここでは内積の実装から、第１引数の長さが優先され、第２引数の配列が第１引数より長い場合でも、 第１引数と同じ長さ分しか積和が計算されない。そのため、配列の長さが不等であることは問題ない。ここで reduceRight の２回めのループが終わり、 akum の内容が [`x_3`] から [`x_3`, `x_2`] となる。

そして３回めの、つまり最後のループを迎える。同様な計算の上、akum の内容として最終的に [`x_3`, `x_2`, `x_1`]が得られる。しかし、 求めているのはこの逆の[`x_1`, `x_2`, `x_3`]なので、最後の reverse() が付け加えられている。

配列もどき

HTML を扱う場合、一見すると配列でありながら、実は配列でない構造がある。 たとえば、getElementsByTagName で取得した DOMエレメントのオブジェクト HTMLCollection や、 querySelectorAll で取得した DOM エレメントのオブジェクト NodeList などの配列風オブジェクト（配列もどき）は、 配列によるアクセスができるので配列と思いがちだが、 実は配列ではない。というのも、配列のクラスに備わっているソートや検索、列挙などが使えないからだ。 

	// link 要素列を取得(HTMLCollection)
	const elements = document.getElementsByTagName('link');
	// 列挙したいがエラー
	elements.forEach(element => console.log(element));

このような場合は、Array.from で配列そのものに変換できる。

	// link 要素列(HTMLCollection)を取得して配列に変換する
	const elements = Array.from(document.getElementsByTagName('link'));
	// 列挙する
	elements.forEach(element => console.log(element));

Array.from() (developer.mozilla.org) 参照。Array.from() は配列もどきの他にも、Set や Map のような反復可能オブジェクトを配列に変換するときにも使える。

別法として、ArrayObject = Array.prototype.slice.call(HTMLCollectionObject) のようにすれば、HTMLCollection が Array に変換できる。 Array.prototype.slice() (developer.mozilla.org) 参照。

数式の記述と表現ほか

本ページで述べた計算は、mathjs で実現されている。 数式記述には ASCIIMathML を使っている。 また数式表現には MathJax4 を使っている。