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

\(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} \) を求める連立一次方程式を解く

連立方程式を解く

  1. 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} \)

  2. 三角行列の行列式は簡単に求めることができるので、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分解

Pivot選択

下三角行列の逆行列

\(\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} \)

上三角行列の逆行列

\(\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} \)

まとめ

なんかできた

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


添付ファイル: filelu2.png 27件 [詳細] filelu.png 28件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS   sitemap
Last-modified: 2017-01-12 (木) 22:37:37 (77d)
ISBN10
ISBN13
9784061426061