こないだのJuliaでNMFしたコードを一般のαダイバージェンスの尺度でNMFできるようにしました． αダイバージェンスの定義はここを参考にしました． 更新式はこの論文の2節に書いてあります．なんと更新式の一部をα倍したり1/α倍したりするだけ，簡単！ using …

各要素が適当な分布に従うテンソルをつくる．例えば，中心がmu, 分散がsigmaの一次元ガウス分布に各要素が従うD×M×Nテンソルをつくりたければ using Distributions mu = 0 sigma = 5 D = 2 M = 3 N = 5 dist = Normal(mu, sigma) X = rand( dist, (D, M, N) )…

集合 B ={ b ∈ {±1}^N } を得たい！ 僕が考えたのはこんな感じ． using Combinatorics function main(N) B=[] for n=1:N combs = collect(combinations(1:N,n)) for comb in combs b = ones(N) for com in comb b[com] = -1 end push!(B,b) end end return B…

自分が勉強するのに役に立ったYoutubeで公開されている講義をどんどん追加していきます． 数学 基礎数学I (1) 多様体論へのいざない 0420 - YouTube 松本多様体の内容を（かなり端折りながらも）丁寧に追っていく授業です．理解が深まりました．面白かった．…

ひとまずpythonにwandbをインストールする． pip install wandb 次に，JuliaにPyCallをインストールする． (@v1.5) pkg> add PyCall wandbの公式ページでアカウントをつくる．適当なプロジェクトを作る． その語，Pythonでwandbを使うようにjuliaで julia > …

julia で string から float への変換は parse() を用いる．単にfloat("0.12")とかはダメ

@time とかいろいろあるみたいですが，私が好きなのはこの方法 using BenchmarkTools runtime = @elapsed begin (処理) end これで変数runtimeに(処理)が終わるまでに要する時間が秒で代入されます． 追記 強い人に教わったのですが result, runtime = @time…

このスライドのp26, アルゴリズム2をJulia言語で実装しました． https://mycourses.aalto.fi/pluginfile.php/1249889/mod_resource/content/1/Lecture6.pdf 素朴なHOSVDや，T-HOSVDより効率よく低ランク近似が実現できているらしいです． ArpackだとフルHOSV…

Julia言語で，T-HOSVDによりテンソルランクを１にするアルゴリズムを実装しました． この論文のアルゴリズム3を実装しました． using TensorToolbox using LinearAlgebra using Arpack function THOSVD(T) # Output is rank-1 tensor N = ndims(T) input_ten…

Julia でSVD をするときに，Lapackを使って using LinearAlgebra n = 10 A = rand(n, n) svd(A) としても良いのですが，これだとO(n3)必要ですね． ランクをrにするのに，完全なSVDをする必要はありません．Arpackのsvdsを使うべし． using Arpack n = 10 r …

pythonのnp.squeezeをjulia言語でやるなら，dropdimsを使うと良い．juliaにも昔は'squeeze'関数があったっぽいが，今は消えてしまった？ docs.julialang.org

collectを使うと良いです． tup = (1,2,3,4,5) ary = collect(tup) これで配列arrayに直されます．

行列の低ランク近似がSVDによって実現されるように，テンソルの低（タッカー）ランク近似がHOSVDで実現できます． なお，Eckart Youngの定理より，SVDがフロベニウスノルムの意味で最良低ランク近似を実現しますが，テンソルの場合はそのような保証はないと…

3×3×3のテンソルを用意したいときは，X = rand(3,3,3) みたいにするとおもいますが，階数が大きい時に X = rand(i1,i2,i3,i4, ..., ,in って手動で入力するのが大変だったので，i = [i1,i2,i3,…,in] を用意して using TensorToolbox X = diagt(i) みたくして…

IBISに参加してみて，「行列分解の研究をしているの行列分解をあんまりよく理解していないな，確率的主成分分析でも勉強するか」と思い，中島変分ベイズの3-5,3-8,6章のほとんど,7章の全てを読んでみた． 須山ベイズを読み込んでいたので，スムーズに読めた…

ただのメモです． 行列Aに対して，exp(A) とすると，行列指数関数が発動されてしまう． element-wiseに演算をしたいときは，.を使う． exp.(A)で各要素ごとにexpされる． A*Bは行列毎の積だけど，A.*Bにすると，要素積になる．

原 啓介「線形性・固有値・テンソル」を読んだ． とてもとても良かった． 物理で線形代数には慣れ親しんでいたはずだけど，いつの間にか私が物理の計算で使うのは線形代数ではなく行列演算になってしまっていて，多重線形性とテンソルの関係とかすっかり忘れ…

楽しい動画だった． 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 =…

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,…

最小値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 =</cmath></iostream></eigen/dense>…

あんまり良い方法が思い浮かばないので，いつもこんな風にやっているのですが，何か良いアイディアがあれば教えてください． #include <Eigen/Dense> #include <iostream> #include <cmath> using namespace Eigen; using namespace std; int matrix_rank(MatrixXd mtx){ FullPivLU<MatrixXd> lu_deco</matrixxd></cmath></iostream></eigen/dense>…

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</fstream></cstdio></time.h></cmath></vector></eigen/dense></iostream>…

先日の記事のコードを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 ; …

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

また強いツイッタラーに教えてもらった． n×n行列のs行目とt列目を削除して，(n-1)×(n-1)行列を得たい，みたいなときに便利なのがInvertedIndices 行列aの取り除きたい列をs，取り除きたい行をtとすると，a[Not(s), Not(t)] で所望の行列を得る． こちらのツ…

ここにある更新式を見て，実装した．白抜きの丸はelement-wiseの積だと理解した． https://www.jjburred.com/research/pdf/jjburred_nmf_updates.pdf収束判定条件はsklearnのNMFの標準設定とできるだけ揃えた． scikit-learn.org以下実装． using LinearAlge…

朝起きて，珈琲飲んで，Surface を見ると電源が切れている（普段はスリープモード）． 変だなと思い，Surfaceを立ち上げてみると，マウス操作のラグがすごい．タッチパネルもすごいラグ．３回再起動したけど治らず．ブルートゥースを切ると症状なくなる．ブ…

機械学習を使ったデジタル・ファブリケーションのためのデザイン支援 （梅谷信行） 空気力学はチャレンジング．空気力学をうまいことやるために，機械学習をつかってみる． 計算空気力学はとても計算コストが高い．計算空気力学は，入力が3D形状で，メッシュ…

JuliaのPlotsとてもカッコ良いのですが，このように軸が見切れたりします． そういうときはusing Plots.PlotMeasuresしたあとに， plot(x,y, bottom_margin = 20mm ) のように，marginオプションを使えば解決できる． left_margin, bottom_margin, right_mar…

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