まぃふぇいばりっと

機械学習やってます.Julia大好きです.勉強したことの殴り書きです.

原 啓介「線形性・固有値・テンソル」は行列代数を卒業して抽象的な線形代数に挑む人におすすめ.

原 啓介「線形性・固有値テンソル」を読んだ.

とてもとても良かった. 物理で線形代数には慣れ親しんでいたはずだけど,いつの間にか私が物理の計算で使うのは線形代数ではなく行列演算になってしまっていて,多重線形性とテンソルの関係とかすっかり忘れてしまっていたので,リハビリとしてもとても良い運動になった..

固有値の計算とか,行列式の計算とか,そういうのに慣れていても,抽象的なベクトル空間の数学としての線形代数の知識が体系的にまとまっていない人に強くおすすめします.テンソルの話や,内積空間の議論も豊富で役に立ちました.

巷の低評価が謎過ぎる!非常に勉強になった.ストーリーがしっかりしていて,何が重要で何が本質なのかが明示的に書かれている.読みやすい.今後も何度も見直すと思う.

階数定理から線形写像の大切な性質が示せるくだりや、行列の積の定義が抽象的な議論から自然に出てくるくだりがとても勉強になった。4章も分かりやすい.線形写像の表現が行列であったように,多重線形写像の表現がテンソル,ふむふむ.(些細なこともしれないけど)双線形写像を線形性の言葉で調べるために,テンソル積を導入するとか,ちゃんと書いてあって本筋を見失わずに読めた.機械学習やってて,急にテンソルを扱うことになり,テンソル積の定義とかをググりながらつまみ食いしがちな私にとっては,とても役に立った本.

同著者の別の本も買ってみようかしら.

Python 極座標上に素数をプロット

楽しい動画だった.

www.youtube.com

練習がてら,任意の素数pに対して極座標(r, θ) = (p, p)の点をプロットするプログラムを作ってみた.

import numpy as np
import matplotlib.pyplot as plt
import math
import time

max_max_number = 100000
interval = 100

def get_prime_numbers():
    start = time.time()
    prime_numbers = [2] 
    for L in range(3, max_max_number):
        chk=True
        for L2 in prime_numbers:
            if L%L2 == 0:chk=False
        if chk==True:prime_numbers.append(L)
    prime_numbers = np.array(prime_numbers)

    return prime_numbers

def converte_pn_coordinates(pns):
    converted_pns_x = [pn*math.cos(pn) for pn in pns]
    converted_pns_y = [pn*math.sin(pn) for pn in pns]
    
    return converted_pns_x, converted_pns_y

def main():
    print("getting prime_numbers...")
    prime_numbers = get_prime_numbers() 
    print("converting polar coordinate...")
    converted_prime_numbers_x, converted_prime_numbers_y = converte_pn_coordinates(prime_numbers)

    for max_number in range(0, max_max_number, interval):
        num_plot_points = len(prime_numbers[prime_numbers < max_number])
        x = converted_prime_numbers_x[:num_plot_points]
        y = converted_prime_numbers_y[:num_plot_points]
        
        max_prime_number = prime_numbers[num_plot_points] 

        plt.figure(figsize=(5,5))
        
        plt.scatter(x, y, s=0.5, c="red")
        plt.title("p < {}".format(max_number))
        plt.axis("off")
        plt.xlim(-max_prime_number, +max_prime_number)
        plt.ylim(-max_prime_number, +max_prime_number)
        plt.savefig("save/{}.png".format(str(max_number).zfill(8)))
        plt.close()

        print("save/{}.png has been saved".format(str(max_number).zfill(8)))

if __name__ == "__main__":[f:id:physics303:20201205224739p:plain]
    main()

学部1年生の時に,リーマン予想とか知らずに,素数の分布を調べようと,1から1000くらいまでの素数をあれこれプロットして規則性がないか調べてたなぁ(何も分からなかった)f:id:physics303:20201205224724p:plain

Python SVDによる低ランク近似

PythonでSVDで入力行列Aのランクをrankに落とすプログラム. ちなみに,C++バージョンはこっち

import numpy as np
import scipy as sp

def SVD(A, rank):
    u, s, v = np.linalg.svd(A)
    ur = u[:, :rank]
    sr = np.matrix(sp.linalg.diagsvd(s[:rank], rank, rank))
    vr = v[:rank, :]

    Rr = ur*sr*vr

    return Rr

ちなみにこのコードだと,O(n3)かかる.Arpack使うと,O(n2 r)でできるはず... scipy.sparse.linalg.svds を参考.(スパース行列じゃなくてもできる)

噂では,sklearn.decomposition.truncatedSVD も使えるようだ. 計算量をさらにO(n2 log r + r2 n)まで落とすこともできるようで,これはこちらの論文に載っている(らしい.私は読んでない)

C++ eigenでランダム行列を作る

最小値LOから最大値HIの値を一様にとる連続分布で,N×Nの行列を作る方法.

#include <Eigen/Dense>
#include <iostream>
#include <cmath>

using namespace Eigen;
using namespace std;

MatrixXd random_uniform_mtx(int N, double LO, double HI){
    double range = HI - LO ;

    MatrixXd P = MatrixXd::Random(N, N);
    P = (P + MatrixXd::Constant(N, N, HI))*range/2;
    P = (P + MatrixXd::Constant(N, N, LO));

    return P;
}

Stack Overflow で見かけたコードだった気がするが,ぐぐっても見当たらなくなっていた...

C++で行列のランクを求める

あんまり良い方法が思い浮かばないので,いつもこんな風にやっているのですが,何か良いアイディアがあれば教えてください.

#include <Eigen/Dense>
#include <iostream>
#include <cmath>

using namespace Eigen;
using namespace std;
int matrix_rank(MatrixXd mtx){
    FullPivLU<MatrixXd> lu_decomp(mtx);
    auto rank = lu_decomp.rank();

    return rank;
}

C++ SVDによる低ランク近似

O(n3)で入力行列Pのランクをrにするアルゴリズムです.eigenをつかっています.

#include <iostream>
#include <Eigen/Dense>
#include <vector>
#include <cmath>
#include <time.h>
#include <cstdio>
#include <fstream>
#include "utils.h"

using namespace std;
using namespace Eigen;

MatrixXd SVD(MatrixXd& P, int r){
    int const N = P.rows();
    int const rank = r;

    //JacobiSVD<MatrixXd> svd(P, ComputeFullU | ComputeFullV);
    BDCSVD<MatrixXd> svd(P, ComputeFullU | ComputeFullV);
    MatrixXd SVD_P(N,N);
    MatrixXd U(N,N), Ur(N,rank-1);
    MatrixXd V(N,N), Vr(N,rank-1);
    MatrixXd S(N,N), Sr(rank-1, rank-1);
    VectorXd sigma(N);

    U = svd.matrixU();
    V = svd.matrixV();
    sigma = svd.singularValues();
    S = sigma.array().matrix().asDiagonal();

    Ur = U(all, seq(0,rank-1));
    Vr = V(all, seq(0,rank-1));
    Sr = S(seq(0,rank-1),seq(0,rank-1));

    SVD_P = Ur * Sr * Vr.transpose();

    return SVD_P;
}

Julia言語 Sinkhorn Knopp アルゴリズムを実装した.

先日の記事のコードをJulia化した.

using LinearAlgebra

function quick_convergence_check(v, v_new, u, u_new, eps)
    if norm( abs.(v - v_new) ) < eps
        if norm( abs.(u - u_new) ) < eps
            return true
        end
    end
    return false
end

function main(A, r, c ; eps = 1.0E-7)
    r = vec(r)
    c = vec(c)
    # vector r and c should be normalized.
    rs = sum(r)
    cs = sum(c)
    r = r / sum(r)
    c = c / sum(c)

    # Initial Step
    v = copy(r)
    u = copy(c)

    cnt = 0
    while true
        # Sinkhorn-knopp Step
        v_new = r ./ ( A' * u )
        u_new = c ./ ( A * v_new )

        if quick_convergence_check(v, v_new, u, u_new, eps)
            break
        end

        # Initial Step
        if cnt > 25000
            println("calculation did not converge")
            exit()
        end

        v = v_new
        u = u_new
        cnt += 1
    end

    P = diagm( vec(u) ) * A * diagm( vec( v ))

    return P
end

A = [1 5 3 6;
     9 5 5 6;
     1 7 10 9;
     2 3 1 4]
r = [7 7 4 2]
c = [1 2 3 4]

P = main(A,r,c)
display(P)
println("\n", sum(P,dims=1) )
println("\n", sum(P,dims=2) )

Python Sinkhorn Knopp アルゴリズムを実装した.

非負行列の行和,列和を揃える操作を行うSinkhorn Knoppアルゴリズムpythonで実装しました.Julia版はこっち. 入力行列Aのi行目の行和がr[i], j列目の列和がc[j]になるようにするアルゴリズムが以下です.

Sinkhorn Knoppアルゴリズムは,エントロピー拘束付き最適化と関りがあります. 最適輸送の文脈での最重要アルゴリズムの1つです.

import numpy as np

def main(A, r, c, eps = 1.0E-7):
    """
    Input  : Matrix A, vectors r and c
             A can be a rectangular matrix
    Output : P = diag(u) A diag(v)
    Sum_i p_1i = c_i / sum(c) 
    Sum_j p_j1 = r_j / sum(r) 
    <=> sum(P.T) = c
        sum(P) = r 
    """

    assert A.ndim == 2, "A is not matrix"
    assert r.ndim == 1 and c.ndim == 1, "r and c should be vectors."
    assert A.shape[0] == c.shape[0], "dim of A does not correspond to r ones."
    assert A.shape[1] == r.shape[0], "dim of A does not correspond to c ones." 

    # vector r and c should be normalized.
    r = r / np.sum(r)
    c = c / np.sum(c)
    
    # Initial Step
    v = np.copy(r) # same as v = r / ( A.T @ c)   
    u = np.copy(c) # same as u = c / ( A @ v)

    cnt = 0
    while True:
        # Sinkhorn-knopp Step
        v_new = r / ( A.T @ u)   
        u_new = c / ( A @ v_new)

        if quick_convergence_check(v, v_new, u, u_new, eps) == True:
            break

        if cnt > 25000:
            print("calculation did not converg")
            exit()
        v = v_new
        u = u_new
        cnt += 1

    P = np.diag(u) @ A @ np.diag(v)
    return P

def quick_convergence_check(v, v_new, u, u_new, eps):
    if np.linalg.norm(abs(v - v_new)) < eps:
        if np.linalg.norm(abs(u - u_new)) < eps:
            return True 
    return False 

if __name__ == "__main__":
    A = np.array([[1,5,3,6],[9,5,5,6],[1,7,10,9],[2,3,1,4]], dtype = "float64")
    c = np.array([7,7,4,2], dtype="float64")
    r = np.array([1,2,3,4], dtype="float64")

    P = main(A, r, c)
    print(sum(P.T))
    print(sum(P))

ちなみに,rとcを全て1からなる行列にすると求まる行列は二重確率行列と呼ばれるのですが,二重確率行列のいろいろな性質についてはこの論文が参考になります.

www.sciencedirect.com

Julia言語 行列の一部を削除する.

また強いツイッタラーに教えてもらった. n×n行列のs行目とt列目を削除して,(n-1)×(n-1)行列を得たい,みたいなときに便利なのがInvertedIndices

行列aの取り除きたい列をs,取り除きたい行をtとすると,a[Not(s), Not(t)] で所望の行列を得る.

f:id:physics303:20201201231949p:plain

こちらのツイートもかなり参考になります.A[begin:end .≠s , begin:end .≠t]でもできる模様

Julia言語 非負値行列因子分解

ここにある更新式を見て,実装した.白抜きの丸はelement-wiseの積だと理解した.
https://www.jjburred.com/research/pdf/jjburred_nmf_updates.pdf

収束判定条件はsklearnのNMFの標準設定とできるだけ揃えた.
scikit-learn.org

以下実装.

using LinearAlgebra
using Random
Random.seed!(123)

function KL(A, B)
    n, m = size(A)
    kl = 0.0
    for i = 1:n
        for j = 1:n
            kl += A[i,j] * log( A[i,j] / B[i,j] ) - A[i,j] + B[i,j]
        end
    end
    return kl
end

function nmf_euc(X, r ; max_iter=200, tol = 0.0001)
    m, n = size(X)
    W = rand(m, r)
    H = rand(r, n)

    error_at_init = norm(X - W*H)
    previous_error = error_at_init
    for iter in 1:max_iter
        H .=  H .* ( W' * X ) ./ ( W' * W * H  )
        W .=  W .* ( X  * H') ./ ( W  * H * H' )

        if tol > 0 && iter % 10 == 0
            error = norm(X - W*H)
            if (previous_error - error) / error_at_init < tol
                break
            end
            previous_error = error
        end
    end

    return W, H
end

function nmf_kl(X, r ; max_iter=200, tol = 0.0001)
    m, n = size(X)
    W = rand(m, r)
    H = rand(r, n)
    one_mn = ones(m, n)

    error_at_init = KL(X, W*H)
    previous_error = error_at_init
    for iter in 1:max_iter
        println(KL(X, W*H))
        H .=  H .* ( W' * ( X ./ (W*H))) ./ ( W' * one_mn )
        W .=  W .* ( (X ./ (W*H) ) * H') ./ ( one_mn * H' )

        if tol > 0 && iter % 10 == 0
            error = KL(X, W*H)
            if (previous_error - error) / error_at_init < tol
                break
            end
            previous_error = error
        end
    end

    return W, H
end

更新するところ,H = hogehoge より,H.= hogehoge の方が微妙に速いって強いツイッタラーに教えてもらった.
Performance Tips · The Julia Language


ちなみに,sklearnでKL-NMFはこんな感じ.デフォルトのsklearnの設定より,上のJulia実装の方が,コスト関数を小さくできる(なぜ?)

import numpy as np
from math import log
from sklearn.decomposition import NMF

rank = 3
nmf = NMF(rank, beta_loss='kullback-leibler', alpha=0, verbose=1,solver='mu') # NMF
W = nmf.fit_transform(X)
H = nmf.components_
X_nmf = np.dot(W, H)


追記

αダイバージェンスでNMFできるように拡張した.
genkaiphd.hatenablog.com

Surface Pro 7 / Bluetoothを ON にすると任意の操作がぐらつく.

朝起きて,珈琲飲んで,Surface を見ると電源が切れている(普段はスリープモード).

 

変だなと思い,Surfaceを立ち上げてみると,マウス操作のラグがすごい.タッチパネルもすごいラグ.3回再起動したけど治らず.ブルートゥースを切ると症状なくなる.ブルートゥースをON/OFF繰り返しても治らない.

 

少し怖いけど,

Bluetoothとその他のデバイス→その他のBluetoothオプション→ハードウェア→Bluetooth Device→プロパティ→ドライバ→デバイスのアンインストール

を実行.次にIntelのサイトからBluetoothドライバをインストールして再起動.

 

これで症状は治った.原因は不明.

 

www.intel.co.jp

IBIS2020 聴講メモ

機械学習を使ったデジタル・ファブリケーションのためのデザイン支援

(梅谷信行)

  • 空気力学はチャレンジング.空気力学をうまいことやるために,機械学習をつかってみる.

  • 計算空気力学はとても計算コストが高い.計算空気力学は,入力が3D形状で,メッシュを生成して,NS方程式をつくって,圧力場や流速場や抵抗係数を求めるが.これが大変.そこで3D形状→回帰分析→抵抗係数にしてしまって,インタラクティブな計算をできるようにした.

  • 訓練データは,実際にいくつかの構造の紙飛行機を飛ばして,その軌跡をトラッキングして得た.

Julia言語 Plotsの軸が見切れる問題

JuliaのPlotsとてもカッコ良いのですが,このように軸が見切れたりします.

f:id:physics303:20201125191441p:plain

そういうときはusing Plots.PlotMeasuresしたあとに, plot(x,y, bottom_margin = 20mm ) のように,marginオプションを使えば解決できる.

left_margin, bottom_margin, right_margin, top_margin など同様に使える.というか単位に"mm"をつけるの,Juliaっぽくていいな.

Julia言語 stringをfloatにする.

x = 3 ; x = float(x) でint型をfloatに変換できますが,x = "3.0" みたいなstring を float(x)としてもうまく変換できません.そういう時は,parse()を使います.

parse(Float64, x)

で無事にstringをFloat64にできます.

Julia言語 コンソールに表を表示させる.

コンソール上にいい感じに数値を並べていくには,以下URLにあるような 

  Using StringLiterals するのがよさげ.  

discourse.julialang.org

ただ,インストールに手間取った. Pkg.clone()はどうも過去の遺物っぽくて,

pkg > add "https://github.com/JuliaString/StringLiterals.jl"

でインストールできた.