Thought GPGPU で逆行列の計算

2019-02-222019
14.4% (52.6/365)
February
77.2% (21.6/28)
Week 08
66.2% (4.6/7)
Day 22 Fri
63.6% (15.2/24)

訪問者

total: 1456
today: 3
yesterday: 4
now: 12

更新

MenuBar

逆行列を求めること＝連立方程式を解くこと
連立方程式を解く
LU分解
Pivot選択
下三角行列の逆行列
上三角行列の逆行列
まとめ
なんかできた

逆行列を求めること＝連立方程式を解くこと †

\(A A^{-1} = E \\ \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{bmatrix} \begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \)

逆行列を求める ⇔ 与えられた \( a_{ij} \) に対する \( x_{ij} \) を求める連立一次方程式を解く

↑

連立方程式を解く †

LU分解する
\(A = LU \)
となる L U を求める
\(\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \\ \end{bmatrix} = \begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix} \begin{bmatrix} 1 & u_{12} & u_{13} & u_{14} \\ 0 & 1 & u_{23} & u_{24} \\ 0 & 0 & 1 & u_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \)
三角行列の行列式は簡単に求めることができるので、L U の逆行列から \(A^{-1}\) を求める
\(A A^{-1} = E \\ L U A^{-1} = E \\ U A^{-1} = L^{-1} E \\ A^{-1} = U^{-1} L^{-1} E //\\ \)

↑

LU分解 †

U の対角成分を 1 に固定して LU を求めてみる
\(\begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix} \begin{bmatrix} 1 & u_{12} & u_{13} & u_{14} \\ 0 & 1 & u_{23} & u_{24} \\ 0 & 0 & 1 & u_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \\ = \begin{bmatrix} l_{11} & l_{11}u_{12} & l_{11}u_{13} & l_{11}u_{14} \\ l_{21} & l_{21}u_{12} + l_{22} & l_{21}u_{13} + l_{22}u_{23} & l_{21}u_{14} + l_{22}u_{24} \\ l_{31} & l_{31}u_{12} + l_{32} & l_{31}u_{13} + l_{32}u_{23} + l_{33} & l_{31}u_{14} + l_{32}u_{24} + l_{33}u_{34} \\ l_{41} & l_{41}u_{12} + l_{42} & l_{41}u_{13} + l_{42}u_{23} + l_{43} & l_{41}u_{14} + l_{42}u_{24} + l_{43}u_{34} + l_{44} \end{bmatrix} \\ = \begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \\ \end{bmatrix} \)
L U は、左上から順に逐次処理で値を決めることができる
\({\color{red}{l_{11}}} = a_{11} \\ l_{11}{\color{red}{u_{12}}} = a_{12} \\ l_{11}{\color{red}{u_{13}}} = a_{13} \\ l_{11}{\color{red}{u_{14}}} = a_{14} \\ {\color{red}{l_{21}}} = a_{21} \\ l_{21}u_{12} + {\color{red}{l_{22}}} = a_{22} \\ l_{21}u_{13} + l_{22}{\color{red}{u_{23}}} = a_{23} \\ l_{21}u_{14} + l_{22}{\color{red}{u_{24}}} = a_{24} \\ {\color{red}{l_{31}}} = a_{31} \\ l_{31}u_{12} + {\color{red}{l_{32}}} = a_{32} \\ l_{31}u_{13} + l_{32}u_{23} + {\color{red}{l_{33}}} = a_{33} \\ l_{31}u_{14} + l_{32}u_{24} + l_{33}{\color{red}{u_{34}}} = a_{34} \\ {\color{red}{l_{41}}} = a_{41} \\ l_{41}u_{12} + {\color{red}{l_{42}}} = a_{42} \\ l_{41}u_{13} + l_{42}u_{23} + {\color{red}{l_{43}}} = a_{43} \\ l_{41}u_{14} + l_{42}u_{24} + l_{43}u_{34} + {\color{red}{l_{44}}} = a_{44} \\ \)
(赤字がその式を評価している時点での未知数 → どの式の未知数も１つなので求めることができる)

\(ie. \\ {\color{red}{l_{ij}}} = a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \\ {\color{red}{u_{ij}}} = \frac{1}{l_{ii}} (a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}) \)
並列度は n-t が限界 (nは列数=行数)
- t行とt列の計算は並列に実行可能。さらに、t行の各要素、t列の各要素は並列に実行可能
- t+1行とt+1列以降の計算は、t行とt列の計算が終わるまで実行できない
- t が変わるごとに並列度が変わるし、外側から玉ねぎを向くように逐次実行が必要 → GPUで計算するのは不得意かな

↑

Pivot選択 †

LU分解では、\(l_{ii}\) が 0 になると \({\color{red}{u_{ij}}}\) の値が求められなくなる
\({\color{red}{l_{ij}}} = a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \\ {\color{red}{u_{ij}}} = \frac{1}{l_{ii}} (a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}) \)
そこで、行列 A の行を適当に入れ替えて、\(l_{ii}\) が 0 にならないようにする。0 にならないようにするばかりではなく、計算誤差を少なくすることを考えると \(|l_{ii}|\) をなるべく大きくしたい。
\(A A^{-1} = E \\ \begin{bmatrix} {\color{red}{a_{11}}} & {\color{red}{a_{12}}} & {\color{red}{a_{13}}} & {\color{red}{a_{14}}} \\ {\color{green}{a_{21}}} & {\color{green}{a_{22}}} & {\color{green}{a_{23}}} & {\color{green}{a_{24}}} \\ {\color{blue}{a_{31}}} & {\color{blue}{a_{32}}} & {\color{blue}{a_{33}}} & {\color{blue}{a_{34}}} \\ {\color{pink}{a_{41}}} & {\color{pink}{a_{42}}} & {\color{pink}{a_{43}}} & {\color{pink}{a_{44}}} \end{bmatrix} \begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{bmatrix} = \begin{bmatrix} {\color{red}{1}} & {\color{red}{0}} & {\color{red}{0}} & {\color{red}{0}} \\ {\color{green}{0}} & {\color{green}{1}} & {\color{green}{0}} & {\color{green}{0}} \\ {\color{blue}{0}} & {\color{blue}{0}} & {\color{blue}{1}} & {\color{blue}{0}} \\ {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{1}} \end{bmatrix} \)

\(A' A^{-1} = B \\ \begin{bmatrix} {\color{blue}{a_{31}}} & {\color{blue}{a_{32}}} & {\color{blue}{a_{33}}} & {\color{blue}{a_{34}}} \\ {\color{green}{a_{21}}} & {\color{green}{a_{22}}} & {\color{green}{a_{23}}} & {\color{green}{a_{24}}} \\ {\color{red}{a_{11}}} & {\color{red}{a_{12}}} & {\color{red}{a_{13}}} & {\color{red}{a_{14}}} \\ {\color{pink}{a_{41}}} & {\color{pink}{a_{42}}} & {\color{pink}{a_{43}}} & {\color{pink}{a_{44}}} \end{bmatrix} \begin{bmatrix} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{bmatrix} = \begin{bmatrix} {\color{blue}{0}} & {\color{blue}{0}} & {\color{blue}{1}} & {\color{blue}{0}} \\ {\color{green}{0}} & {\color{green}{1}} & {\color{green}{0}} & {\color{green}{0}} \\ {\color{red}{1}} & {\color{red}{0}} & {\color{red}{0}} & {\color{red}{0}} \\ {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{1}} \end{bmatrix} \)
LU分解の式を修正
\(A A^{-1} = E \\ A' A^{-1} = B \\ L' U' A^{-1} = B \\ U' A^{-1} = L'^{-1} B \\ A^{-1} = U'^{-1} L'^{-1} B \\ \)
一般の一次方程式を解くときには、B を掛ける必要があるけど、逆行列を求めるとき限定で \(Y = U'^{-1} L'^{-1}\) の列を入れ替えてやるだけでいい
\(A^{-1} = U'^{-1} L'^{-1} B = Y B = \\ \begin{bmatrix} y_{11} & y_{12} & y_{13} & y_{14} \\ y_{21} & y_{22} & y_{23} & y_{24} \\ y_{31} & y_{32} & y_{33} & y_{34} \\ y_{41} & y_{42} & y_{43} & y_{44} \end{bmatrix} \begin{bmatrix} {\color{blue}{0}} & {\color{blue}{0}} & {\color{blue}{1}} & {\color{blue}{0}} \\ {\color{green}{0}} & {\color{green}{1}} & {\color{green}{0}} & {\color{green}{0}} \\ {\color{red}{1}} & {\color{red}{0}} & {\color{red}{0}} & {\color{red}{0}} \\ {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{0}} & {\color{pink}{1}} \end{bmatrix} = \begin{bmatrix} {\color{red}{y_{13}}} & {\color{green}{y_{12}}} & {\color{blue}{y_{11}}} & {\color{pink}{y_{14}}} \\ {\color{red}{y_{23}}} & {\color{green}{y_{22}}} & {\color{blue}{y_{21}}} & {\color{pink}{y_{24}}} \\ {\color{red}{y_{33}}} & {\color{green}{y_{32}}} & {\color{blue}{y_{31}}} & {\color{pink}{y_{34}}} \\ {\color{red}{y_{43}}} & {\color{green}{y_{42}}} & {\color{blue}{y_{41}}} & {\color{pink}{y_{44}}} \end{bmatrix} \)
「行列 A の行を適当に入れ替えて」ってどうすんの？
1. 未計算分の \(l_{ij}\) を縦方向に計算する
  \({\color{red}{l_{ij}}} = a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj} \\ \)
2. 一番大きい \(l_{ij}\) を持つ行を未計算部分の一番上の行に入れ替える (=Pivot選択)
3. 選んだ Pivot 行の \(l_{ij}\) を横方向に計算する
  \({\color{red}{u_{ij}}} = \frac{1}{l_{ii}} (a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}) \)
並列に実行できる計算
- 上図の③と④の部分 (③、④の各要素の計算は並列実行可能)
- Pivot行の \(u_{t\ t+1}\) を計算したあとの
  - \(u_{t\ j}\) j=t+2,t+3,...,n の計算
  - \(u_{i\ t+1}\) i=t+1,t+2,...,n の計算と次の Pivot 行の選択
- CPU で計算するなら、④をメインスレッドにして、③をサブスレッドにする。④(=①+②) の処理終了後に③の終了待ち(JOIN)を入れることになるかな
ピボット選択をしても \(l_{ii}=0\) となる場合は、逆行列を計算できない行列だった
今回やった行方向のピボット選択に加えて、列方向のピボット選択もできるが、処理がややこしいことになるんで普通はやらない
逆行列の計算では、巨大な行列を扱うことが有るので、とにかくワークエリア無しですむように実装する必要あり

↑

下三角行列の逆行列 †

\(\begin{bmatrix} l_{11} & 0 & 0 & 0 \\ l_{21} & l_{22} & 0 & 0 \\ l_{31} & l_{32} & l_{33} & 0 \\ l_{41} & l_{42} & l_{43} & l_{44} \end{bmatrix} \begin{bmatrix} l^{-1}_{11} & 0 & 0 & 0 \\ l^{-1}_{21} & l^{-1}_{22} & 0 & 0 \\ l^{-1}_{31} & l^{-1}_{32} & l^{-1}_{33} & 0 \\ l^{-1}_{41} & l^{-1}_{42} & l^{-1}_{43} & l^{-1}_{44} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \)

1列目
\(l_{11} l^{-1}_{11} = 1 \\ l_{21} l^{-1}_{11} + l_{22} l^{-1}_{21} = 0 \\ l_{31} l^{-1}_{11} + l_{32} l^{-1}_{21} + l_{33} l^{-1}_{31} = 0 \\ l_{41} l^{-1}_{11} + l_{42} l^{-1}_{21} + l_{43} l^{-1}_{31} + l_{44} l^{-1}_{41} = 0 \)
ie.
\(l^{-1}_{11} = \frac{1}{l_{11}} \\ l^{-1}_{21} = -\frac{1}{l_{22}}(l_{21} l^{-1}_{11}) \\ l^{-1}_{31} = -\frac{1}{l_{33}}(l_{31} l^{-1}_{11} + l_{32} l^{-1}_{21}) \\ l^{-1}_{41} = -\frac{1}{l_{44}}(l_{41} l^{-1}_{11} + l_{42} l^{-1}_{21} + l_{43} l^{-1}_{31}) \)
2列目
\(l_{22} l^{-1}_{22} = 1 \\ l_{32} l^{-1}_{22} + l_{33} l^{-1}_{32} = 0 \\ l_{42} l^{-1}_{22} + l_{43} l^{-1}_{32} + l_{44} l^{-1}_{42} = 0 \)
ie.
\(l^{-1}_{22} = \frac{1}{l_{22}} \\ l^{-1}_{32} = -\frac{1}{l_{33}}(l_{32} l^{-1}_{22}) \\ l^{-1}_{42} = -\frac{1}{l_{44}}(l_{42} l^{-1}_{22} + l_{43} l^{-1}_{32}) \)
3列目
\(l_{33} l^{-1}_{33} = 1 \\ l_{43} l^{-1}_{33} + l_{44} l^{-1}_{43} = 0 \)
ie.
\(l^{-1}_{33} = \frac{1}{l_{33}} \\ l^{-1}_{43} = -\frac{1}{l_{44}}(l_{43} l^{-1}_{33}) \)
4列目
\(l_{44} l^{-1}_{44} = 1 \)
ie.
\(l^{-1}_{44} = \frac{1}{l_{44}} \)
一般化 k 列目
\(l^{-1}_{kk} = \frac{1}{l_{kk}} \\ l^{-1}_{ik} = -\frac{1}{l_{ii}} \sum_{j=k}^{i-1} l_{ij} l^{-1}_{jk} \\ (i=k+1,k+2,...,n) \)
並列度
- 元データ格納領域とは別に、逆行列の格納領域を作るのであれば、各列並列で実行できる → 並列度 n (行数)
- 元データ格納領域に、逆行列を格納する場合には並列実行できない（逐次処理になる)

↑

上三角行列の逆行列 †

\(\begin{bmatrix} 1 & u_{12} & u_{13} & u_{14} \\ 0 & 1 & u_{23} & u_{24} \\ 0 & 0 & 1 & u_{34} \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} u^{-1}_{11} & u^{-1}_{12} & u^{-1}_{13} & u^{-1}_{14} \\ 0 & u^{-1}_{22} & u^{-1}_{23} & u^{-1}_{24} \\ 0 & 0 & u^{-1}_{33} & u^{-1}_{34} \\ 0 & 0 & 0 & u^{-1}_{44} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \)

1列目
\(u^{-1}_{11} = 1 \)
2列目
\(u^{-1}_{12} + u_{12}u^{-1}_{22} = 0 \\ u^{-1}_{22} = 1 \)
ie.
\(u^{-1}_{22} = 1 \\ u^{-1}_{12} = -(u_{12} u^{-1}_{22}) \)
3列目
\(u^{-1}_{13} + u_{12} u^{-1}_{23} + u_{13}u^{-1}_{33} = 0 \\ u^{-1}_{23} + u_{23} u^{-1}_{33} = 0 \\ u^{-1}_{33} = 1 \)
ie.
\(u^{-1}_{33} = 1 \\ u^{-1}_{23} = -(u_{23} u^{-1}_{33}) \\ u^{-1}_{13} = -(u_{12} u^{-1}_{23} + u_{13}u^{-1}_{33}) \)
4列目
\(u^{-1}_{14} + u_{12} u^{-1}_{24} + u_{13} u^{-1}_{34} + u_{14}u^{-1}_{44} = 0 \\ u^{-1}_{24} + u_{23} u^{-1}_{34} + u_{24}u^{-1}_{44} = 0 \\ u^{-1}_{34} + u_{34}u^{-1}_{44} = 0 \\ u^{-1}_{44} = 1 \)
ie.
\(u^{-1}_{44} = 1 \\ u^{-1}_{34} = -(u_{34} u^{-1}_{44}) \\ u^{-1}_{24} = -(u_{23} u^{-1}_{34} + u_{24} u^{-1}_{44}) \\ u^{-1}_{14} = -(u_{12} u^{-1}_{24} + u_{13} u^{-1}_{34} + u_{14} u^{-1}_{44}) \)
一般化 k 列目
\(u^{-1}_{kk}=1 \\ u^{-1}_{ik} = -\sum_{j=i+1}^{k} u_{ij} u^{-1}_{jk} \\ (i=k-1,k-2,...,1) \)
並列度
- 元データ格納領域とは別に、逆行列の格納領域を作るのであれば、各列並列で実行できる → 並列度 n (行数)
- 下三角行列の列を求める計算の続きで計算すると、各列の計算量がおなじになる
- 元データ格納領域に、逆行列を格納する場合には並列実行できない（逐次処理になる)

↑

まとめ †

LU分解の並列度は n-k (n:列数=行数, k=1,2,3,...,n)。CPUでやったほうが良いかな
LとUの逆行列を求める計算の並列度は n (列数)。CPU でやるか GPU でやるかは微妙なところ
U^-1L^-1 の計算は GPU の出番

↑

なんかできた †

https://github.com/kagyuu/Matarapi/blob/master/src/main/java/com/matarapi/MatKernel.java

/*
 * Copyright 2017 HONDOH Atsushi.
 */
package com.matarapi;

import com.aparapi.Kernel;

/**
 * Matrix calculate Kernel.
 * @author a.ho
 */
public abstract class MatKernel extends Kernel {
    /**
     * data array.
     */
    protected float ary[];
    /**
     * data offset.
     */
    protected int[] offset;
    /**
     * matrix size.
     */
    protected int[] matSize;
    /**
     * col size.
     */
    protected int[] colSize;

    /**
     * Constructor.
     * @param data matrixes
     */
    public MatKernel(IMatrix ...data) {
        int arySize = 0;
        for (IMatrix m : data) {
            arySize += m.getSize();
        }
        ary = new float[arySize];
        offset = new int[data.length];
        matSize = new int[data.length];
        colSize = new int[data.length];

        int p = 0;
        int n = 0;
        for (IMatrix m : data) {

            offset[n] = p;
            matSize[n] = m.getRowSize() * m.getColSize();
            colSize[n] = m.getColSize();

            p = m.accept(ary, p);
            n += 1;
        }
    }

    /**
     * get Matrix.
     * @param no matrix number
     * @return Matrix
     */
    public Matrix getMat(final int no) {
        return new Matrix(ary, offset[no], matSize[no] / colSize[no], colSize[no]);
    }

    /**
     * add.
     * <pre>
     * in1       in2        out
     * |a b c|   |A B C|    |a+A b+B c+C|
     * |d e f| + |D E F| -> |d+D e+E f+F|
     * |g h i|   |G H I|    |g+G h+H i+I|
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in1 input1
     * @param in2 input2
     * @param out output
     */
    protected void matAdd(int in1, int in2, int out) {

        int i = getGlobalId();
        if (i > matSize[out]) return;

        int p0 = offset[out] + i;
        int p1 = offset[in1] + i;
        int p2 = offset[in2] + i;
        ary[p0] = ary[p1] + ary[p2];
    }

    /**
     * sub.
     * <pre>
     * in1       in2        out
     * |a b c|   |A B C|    |a-A b-B c-C|
     * |d e f| - |D E F| -> |d-D e-E f-F|
     * |g h i|   |G H I|    |g-G h-H i-I|
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in1 input1
     * @param in2 input2
     * @param out output
     */
    protected void matSub(int in1, int in2, int out) {

        int i = getGlobalId();
        if (i > matSize[out]) return;

        int p0 = offset[out] + i;
        int p1 = offset[in1] + i;
        int p2 = offset[in2] + i;
        ary[p0] = ary[p1] - ary[p2];
    }

    /**
     * ReLU (Rectified linear unit).
     * <pre>
     * An simple activation function for neural network.
     * This is enough complicity for machine learning.
     * That's differentiate is so simple!
     *
     * ReLu(x) = max(0,x)
     *         = 0 (x <  0)
     *           x (x >= 0)
     *
     * cf.
     * ReLu-1(x) = 0 (x <  0)
     *             1 (x >= 0)
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in input
     * @param out output
     */
    protected void matReLu(int in,int out) {

        int i = getGlobalId();
        if (i > matSize[out]) return;

        int p0 = offset[out] + i;
        int p1 = offset[in] + i;
        ary[p0] = max(ary[p1],0.0f);
    }

    /**
     * differentiate of ReLU (Rectified linear unit).
     * <pre>
     * an simple activation function for neural network.
     *
     * ReLu-1(x) = 0 (x <  0)
     *             1 (x >= 0)
     * 
     * cf.
     * ReLu(x) = max(0,x)
     *         = 0 (x <  0)
     *           x (x >= 0)
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in input
     * @param out output
     */
    protected void matdReLu(int in,int out) {
        
        int i = getGlobalId();
        if (i > matSize[out]) return;

        int p0 = offset[out] + i;
        int p1 = offset[in] + i;
        ary[p0] = ary[p1] < 0.0f ? 0.0f : 1.0f;
    }

    /**
     * multiple.
     * <pre>
     * in1       in2        out
     * |a b c|   |A B C|    |aA+bD+cG dA+eD+fG gA+hD+iG|
     * |d e f| X |D E F| -> |aB+bE+cH dB+eE+fH gB+hE+iH|
     * |g h i|   |G H I|    |aC+aF+aI dC+eF+fI gC+hF+iI|
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in1 input1
     * @param in2 input2
     * @param out output
     */
    protected void matMul(int in1, int in2, int out) {
        
        int i = getGlobalId();
        if (i > matSize[out]) return;

        int c0 = colSize[out];
        int c1 = colSize[in1];
        int c2 = colSize[in2];

        int col = i % c0;
        int row = (i - col) / c0;

        int p0 = offset[out] + i;
        int p1 = offset[in1] + c1 * row;
        int p2 = offset[in2] + col;

        float mul = 0.0f;
        for (int cnt=0; cnt < c1; cnt++) {
            mul += ary[p1] * ary[p2];
            p1 += 1;
            p2 += c2;
        }
        ary[p0] = mul;
    }
    
    /**
     * multiple LU.
     * <pre>
     * in1       in2        out
     * |a 0 0|   |1 B C|    |a aB   aC      |
     * |d e 0| X |0 1 F| -> |d dB+e dC+eF   |
     * |g h i|   |0 0 1|    |g gB+h gC+hF+i |
     * 
     * in1 and in2 is combined represnted in arg 'in'
     * |a B C|
     * |d e F|
     * |g h i|
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in  input
     * @param out output
     */
    protected void matMulLU(int in, int out) {
        
        int i = getGlobalId();
        if (i > matSize[out]) return;
        
        int size = colSize[in];

        int col = i % size;
        int row = (i - col) / size;
        
        int p0 = offset[out] + i;
        int pL = offset[in] + size * row;
        int pU = offset[in] + col;

        float mul = 0.0f;
        int step = min(col, row);
        for (int cnt=0; cnt <= step; cnt++) {
            if (col == cnt) {
                mul += ary[pL];
            } else {
                mul += ary[pL] * ary[pU];
            }
            pL += 1;
            pU += size;
        }
        ary[p0] = mul;
    }

    /**
     * multiple UL.
     * <pre>
     * in1       in2        out
     * |1 B C|   |a 0 0|    |a+Bd+Cg Be+Ch Ci |
     * |0 1 F| X |d e 0| -> |d+Fg    e+FH  Fi |
     * |0 0 1|   |g h i|    |g       h     i  |
     * 
     * in1 and in2 is combined represnted in arg 'in'
     * |a B C|
     * |d e F|
     * |g h i|
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in  input
     * @param out output
     */
    protected void matMulUL(int in, int out) {
        
        int i = getGlobalId();
        if (i > matSize[out]) return;
        
        int size = colSize[in];

        int col = i % size;
        int row = (i - col) / size;
        int skip = max(row, col);
        
        int p0 = offset[out] + i;
        int pU = offset[in] + size * row + skip;
        int pL = offset[in] + col + size * skip;

        float mul = 0.0f;
        for (int cnt=skip; cnt < size; cnt++) {
            if (row == cnt) {
                mul += ary[pL];
            } else {
                mul += ary[pU] * ary[pL];
            }
            pU += 1;
            pL += size;
        }
        ary[p0] = mul;
    }
    
    /**
     * ivert Upper triangle matrix.
     * <pre>
     * in        out        E
     * |1 b c|   |1 B C|    |1 0 0|
     * |0 1 f| X |0 1 F| -> |0 1 0|
     * |0 0 1|   |0 0 1|    |0 0 1|
     * 
     * The Lower triangle of input matrix, includes diagonal components '1'
     * , is ignored.
     * 
     * paralles size : col size
     * </pre>
     * @param in  input
     * @param out output
     */
    protected void matInvU(int in, int out) {
        
        int k = getGlobalId();
        
        int size = colSize[in];
        
        if (k > size) { return; }
        
        int p0 = offset[out] + size * k + k;        
        ary[p0] = 1.0f;
        
        float inv;
        int pU, pUi;
        for (int i = k - 1; i >= 0; i--) {
            inv = 0.0f;
            
            int j = i + 1;
            pU = offset[in] + size * i + j; // U[i][j]
            pUi = offset[out] + size * j + k; // U-1[j][k]
            for (; j <= k; j++) {
                inv -= ary[pU] * ary[pUi];
                
                pU += 1;     // move right
                pUi += size; // move down
            }
            p0 -= size; // move up
            ary[p0] = inv;
        }
    }
    
    /**
     * ivert Lower triangle matrix.
     * <pre>
     * in        out        E
     * |a 0 0|   |A 0 0|    |1 0 0|
     * |d e 0| X |D E 0| -> |0 1 0|
     * |g h i|   |G H I|    |0 0 1|
     * 
     * The Upper triangle of input matrix is ignored.
     * 
     * paralles size : col size
     * </pre>
     * @param in  input
     * @param out output
     */
    protected void matInvL(int in, int out) {
        
        int k = getGlobalId();
        
        int size = colSize[in];
        
        if (k > size) { return; }
        
        int p0 = offset[out] + size * k + k;
        int p1 = offset[in] + size * k + k;
        ary[p0] = 1.0f / ary[p1];
        
        float inv;
        int pL, pLi;
        for (int i = k + 1; i < size; i++) {
            inv = 0.0f;
            
            int j = k;
            pL = offset[in] + size * i + j; // L[i][j]
            pLi = offset[out] + size * j + j; // L-1[j][k]
            for (; j <= i-1; j++) {
                inv -= ary[pL] * ary[pLi];
                
                pL += 1;     // move right
                pLi += size; // move down
            }
            p0 += size; // move down
            // L-1[i][k] = - (Σ(L[i][j]*L-1[j][k]) / L[i][i]
            ary[p0] = inv / ary[offset[in] + size * i + i]; 
        }
    }
    
    /**
     * calculate hadamard product.
     * <pre>
     * in1       in2        out
     * |a b c|   |A B C|    |aA bB cC|
     * |d e f| O |D E F| -> |dD eE fF|
     * |g h i|   |G H I|    |gG hH iI|
     * </pre>
     * @param in1 input1
     * @param in2 input2
     * @param out output
     */
    protected void matHmul(int in1, int in2, int out) {

        int i = getGlobalId();
        if (i > matSize[out]) return;

        int p0 = offset[out] + i;
        int p1 = offset[in1] + i;
        int p2 = offset[in2] + i;
        ary[p0] = ary[p1] * ary[p2];
    }

    /**
     * calculate transpose matrix.
     * <pre>
     * in        out
     * |a b c| = |a d|
     * |d e f|   |b e|
     *           |c f|
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in input
     * @param out output
     */
    protected void matTranspose(int in, int out) {

        int i = getGlobalId();
        if (i > matSize[out]) return;

        int c0 = colSize[out];
        int c1 = colSize[in];

        int col = i % c1;
        int row = (i - col) / c1;

        int p1 = offset[in] + i;
        int p0 = offset[out] + row + col * c0; // col <-> row

        ary[p0] = ary[p1];
    }

    /**
     * sort columns.
     * <pre>
     * in        out
     * |a b c| = |b a c|
     * |d e f|   |e d f|
     * |g h i|   |h g i|
     * 
     * order = [1,0,2] means:
     * Col 0 -> Col 1
     * Col 1 -> Col 0
     * Col 2 -> Col 2
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in input
     * @param out output
     * @param order sort order
     */
    protected void matSortCol(int in, int out, int[] order) {

        int i = getGlobalId();
        if (i > matSize[out]) return;

        int c0 = colSize[out];
        int c1 = colSize[in];

        int col = i % c1;
        int row = (i - col) / c1;
        
        int p0 = offset[out] + order[col] + c0 * row;
        int p1 = offset[in] + i;

        ary[p0] = ary[p1];
    }
    
    /**
     * sort rows.
     * <pre>
     * in        out
     * |a b c| = |d e f|
     * |d e f|   |a b c|
     * |g h i|   |g h i|
     * 
     * order = [1,0,2] means:
     * Row 0 -> Row 1
     * Row 1 -> Row 0
     * Row 2 -> Row 2
     * 
     * paralles size : matrix size (= row * col)
     * </pre>
     * @param in input
     * @param out output
     * @param order sort order
     */
    protected void matSortRow(int in, int out, int[] order) {

        int i = getGlobalId();
        if (i > matSize[out]) return;

        int c0 = colSize[out];
        int c1 = colSize[in];

        int col = i % c1;
        int row = (i - col) / c1;

        int p0 = offset[out] + col + c0 * order[row];
        int p1 = offset[in] + i;

        ary[p0] = ary[p1];
    }

    @Override
    public String toString() {
        StringBuilder sb = new StringBuilder();

        sb.append("ary:\n");
        for (float val : ary) {
            sb.append(String.format("%.2f ", val));
        }
        sb.append("\noffset:\n");
        for (int val : offset) {
            sb.append(String.format("%d ", val));
        }
        sb.append("\nrow size:\n");
        int cnt = 0;
        for (int val : matSize) {
            sb.append(String.format("%d ", val / colSize[cnt++]));
        }
        sb.append("\ncol size:\n");
        for (int val : colSize) {
            sb.append(String.format("%d ", val));
        }

        return sb.toString();
    }
}

https://github.com/kagyuu/Matarapi/blob/master/src/test/java/com/matarapi/test/MatrixTest.java

/*
 * Copyright 2017 HONDOH Atsushi.
 */
package com.matarapi.test;

import com.matarapi.IMatrix;
import com.matarapi.MatKernel;
import com.matarapi.Matrix;
import com.matarapi.MockMatrix;
import static com.matarapi.test.MatAssert.assertMatrix;
import org.junit.Test;
import static com.matarapi.test.MatAssert.assertMatrixInvert;

/**
 * Simple calculation test.
 *
 * @author atsushi
 */
public class MatrixTest {

    @Test
    public void testLU() throws Exception {
        Matrix m0 = new Matrix(new float[][]{
            {2.0f, 5.0f, 7.0f, 8.0f},
            {4.0f, 13.0f, 20.0f, 25.0f},
            {8.0f, 29.0f, 50.0f, 71.0f},
            {10.0f, 34.0f, 78.0f, 98.0f},
        });
        // m0 will be over-writed to lu-matrix
        Matrix original = m0.copyAll();
        
        // m0 => LU
        // oreder is row order for LU dividing
        final int[] order = m0.toLU();
                        
        // L U => mul
        IMatrix work = new MockMatrix(4, 4);
        MatKernel kernel = new MatKernel(m0, work) {

            @Override
            public void run() {
                int pass = getPassId();
                if (0 == pass) {
                    matMulLU(0, 1);
                } else if (1 == pass) {
                    // revert pivoting
                    matSortRow(1, 0, order);
                }
            }
        };
        kernel.execute(m0.getSize(),2);
        Matrix mul = kernel.getMat(0);
                
//        System.out.println(mul.toString());
        
        // LU expects m0
        assertMatrix(original, mul, 1.0E-6f);        
    }
    
    @Test
    public void testInvert() throws Exception {
        Matrix m0 = new Matrix(new float[][]{
            {2.0f, 5.0f, 7.0f, 8.0f},
            {4.0f, 13.0f, 20.0f, 25.0f},
            {8.0f, 29.0f, 50.0f, 71.0f},
            {10.0f, 34.0f, 78.0f, 98.0f},
        });
        // m0 will be over-writed to lu-matrix
        Matrix original = m0.copyAll();
        
        // m0 => LU with pivoting
        // oreder is row order for LU dividing
        final int[] order = m0.toLU();
                                
        // L U => mul
        IMatrix work = new MockMatrix(4, 4);
        MatKernel kernel = new MatKernel(m0, work) {

            @Override
            public void run() {
                int pass = getPassId();
                // Can't use switch statement on Aparapi.
                if (0 == pass) {
                    // U => U(-1)
                    matInvU(0,1);
                    // L => L(-1)
                    matInvL(0,1);
                } else if (1 == pass) {
                    // U(-1) L(-1) => A(-1)
                    matMulUL(1,0);
                } else if (2 == pass) {
                    // revert pivoting                    
                    matSortCol(0,1,order);
                }
            }
        };
        kernel.execute(m0.getSize(), 3);
        Matrix invert = kernel.getMat(1);
        
//        System.out.println(expect);
//        System.out.println(invert);
        
        // precision is not well, but its ok.
        // It's enough for the machine learning.
        // We don't use GPGPU for Rocket Science.
        assertMatrixInvert(original, invert, 1.0E-5f);
    }
}

Culture

Thought GPGPU で逆行列の計算

目次

訪問者

更新

最新の10件

人気の10件

MenuBar

逆行列を求めること＝連立方程式を解くこと †

連立方程式を解く †

LU分解 †

Pivot選択 †

下三角行列の逆行列 †

上三角行列の逆行列 †

まとめ †

なんかできた †

https://github.com/kagyuu/Matarapi/blob/master/src/main/java/com/matarapi/MatKernel.java

https://github.com/kagyuu/Matarapi/blob/master/src/test/java/com/matarapi/test/MatrixTest.java

Thought GPGPU で逆行列の計算

目次

訪問者

更新

最新の10件

人気の10件

MenuBar

逆行列を求めること ＝ 連立方程式を解くこと †

連立方程式を解く †

LU分解 †

Pivot選択 †

下三角行列の逆行列 †

上三角行列の逆行列 †

まとめ †

なんかできた †

https://github.com/kagyuu/Matarapi/blob/master/src/main/java/com/matarapi/MatKernel.java

https://github.com/kagyuu/Matarapi/blob/master/src/test/java/com/matarapi/test/MatrixTest.java

逆行列を求めること＝連立方程式を解くこと †