KLダイバージェンスを最小化するNMFはよく知られていますし，sklearnで'beta_loss=1'にすればすぐに利用できます． 今回はKLダイバージェンスの拡張であるαダイバージェンスを小さくするNMFを実装しました． この論文の更新式をさくっと実装しました． www.s…

多変数関数f(x)を確率単体上のxで最適化したい．（最近，情報幾何をやっているので，こういうタスクばかりやってる気がする） つまり，xは 0 < x[i] < 1 かつ sum(x) = 1 の範囲でf(x)の最小化をしたい． ちまちま自分でソルバー書いてたけど，scipy使う方が…

BenchmarkToolsって便利ですよね．例のように(10×10)のランダム行列の固有値を計算するのに要した時間を計測したとします． ここから，最小時間37.748とか，中央値47.454とかを取り出すにはどうしたらよいでしょうか，というのがいつも分からなくなるのでメ…

しょっちゅう分からなくなってはTwitterで強い人に教えてもらうので，ここにメモ． CartesianIndexは.Iで添え字を取り出せる．これを'collect'すると所望のvecになる. ただ，これだとちょっと気持ち悪いと思う人もいるようなので， Tuple(CartesianIndex(2,3…

前回，補助行列が一つだけのNMMFの実装を公開した． 今回は同様に，これらの論文にのっとり，補助行列が２つある場合のNMMFを実装した． 行列X,Y,Zを，WH, AH,WBに同時に分解するというタスクである．コスト関数はKL情報量．（LSの場合の更新式ってどこかに…

複合データ分析技術についてこのpdfで勉強してたら，NMMFというものを知った． ちゃんとした議論は，IJCAIでしている． www.ijcai.org 人工知能学会の予稿や論文，絵本のパターン抽出や，和文論文が参考になる． 私も研究で使うことになったのでささっと実装…

昨日の記事の完全版．Nonnegative Tucker Decomposition の Julia 実装． 重要なアルゴリズムと思うけど，Julia版は探しても見つからなかったので，自分で書いた． ieeexplore.ieee.org Python版はすでに実装があるみたい．これはSVDで初期化しているっぽい…

↓ こっちのほうが役に立つと思う genkaiphd.hatenablog.com ↑ こっちのほうが役に立つと思う 研究に必要なので，非負テンソルのタッカー分解の論文をjuliaで実装した．とりあえず，コスト関数がLSエラーのものを実装した．KLダイバージェンスの場合について…

一番手っ取り早いのはこんな感じ A = rand(5,5) A[A.<1.0e-10] .= 0 以下，強い人に教わったこと． 絶対値が1e-10以下なら0にする場合に真っ先に思い浮かぶのは X[abs.(X) .≤ 1e-10] .= 0 だけど，これは遅い． やりたいこと分解して， small2zero(x, ε=1e-1…

重み付き非負行列分解というタスクがある．欠損ありのNMFとかに使える． 研究に必要だったので，オーソドックスな次の手法を実装した． www.semanticscholar.org 実装は簡単で普通のNMFの更新式をちょっと変えるだけ． フロベニウスノルム尺度のものとKL尺度…

左確率行列分解（Left-Stochastic Matrix Factorization）というタスクがある．NMFに対称性と列規格化を要請する．これでうまいことクラスタリングができる． 研究で比較実験を行うためにクラスタ数2の場合のみを実装した． jmlr.org アルゴリズム中に確率単…

修士課程までディラック電子系の理論物理をやってました．そのまま博士課程に進学することも考えましたが，色々考えて民間に就職して，深層学習を使ったコンピュタービジョンの研究をしました．２年後，会社を辞めて情報系の博士課程に編入しました．先日，…

非負の対称行列を低ランク二重確率行列で近似するタスクに関する論文を読んだ． ここでいう近似は，KL情報量の意味での近似． jmlr.org 論文中のAlgorithm 1を実装した． using LinearAlgebra function D_KL(A,B) (N, M) = size(A) dkl = 0.0 for i=1:N for …

emps = [ [], [], [], [], [] ] みたいなのをつくりたい．そういう時は， N = 5 emps = [[] for i=1:N] とすればよい． An array of empty arrays - First steps - JuliaLang

大昔に研究で使っていたCのコードが発掘された． グラフェンのジグザグナノリボンのバンド図を求めます。Nはユニットセル内の原子の数です。最近接のとびうつりのみを考えたタイトバインディングで計算しています。計算結果は/dataにdatファイルとして出力し…

六方格子系に垂直に磁場を与えることでHofstadter's butterflyが現れる． 昔かいたFortranコードを発掘したのでアップロード． 自然界に現れるフラクタル構造って良いですよね． ちなみに，Hofstadter's butterflyは六方格子系以外の系に垂直磁場を与えても…

すっごい勉強になった．自分の専門分野もろかぶりだった． 情報幾何と機械学習あたりやっている人にとっては役に立つと思う． 特に印象に残ってるのは次の３つ． 江口真透：情報幾何学入門 1.通常の最小二乗法はユークリッド幾何の三平方の定理 2.目的変数に…

これを見つけるのに時間がかかってしまった． \usepackage{xr}を使えば良いらしいというのは，このページをみればわかったのだが，overleafでうまくいかなかった．いろいろググると，次のような解説を見つけた．これでいけた！どうやらなぞのlatexmkrcを作る…

研究に必要だったので，Julia言語でこの論文のアルゴリズム3を実装した．「KL情報量でCP分解しましょう」って話だけど，中身はポアソン分布とかいろいろ書いてあってちゃんと読んでないから分からん．とりあえず非負のテンソルCP分解の手法の一つ． arxiv.or…

fMRIの出力がniiファイルらしい．良く分からないが，using NIfTI で開ける． using NIfTI ni = niread("hogehoge.nii") ni = niread("foo.nii.gz") でOK．普通の配列にしたかったら convert(Array{Float64}, ni) でOK． github.com

NASAのページからgeotifファイルをダウンロードできるんだけど，juliaでどう開けばよいか分からなかった． 色々調べたけど，多分これが一番良い． github.com using GeoArrays geoarray = GeoArrays.read("hogehoge.tif") 座標をとりだしたかったら， coords…

juliaで using FileIO img = load("hogehoge.png") したあとに，r,g,b要素を取り出したかったら，それぞれ using Colors r = red.(img) g = green.(img) b = blue.(img) とする．ちなみにこのままだとr,g,bはそれぞれ Array{N0f8,2} with eltype FixedPointN…

自然数nに対して，[n] = {1, 2, …, n} として，[I]×[J]×[K]が欲しい時ととかってあるじゃないですか．テンソルの添え字でloop作りたいときとか．そういうときは，IterToolsのproductを使うとよいことに気が付いた． using IterTools I = 2 J = 3 K = 2 for i…

社会人学生 Advent Calendar 2022 博士論文執筆の際にお願いしたいこと 業績ゼロから始める社会人博士 理工系大学・高専の研究室不登校 (修士/博士/普通の)論文執筆の際にお願いしたいこと、その二 A guide to Scientifc Research 博士課程後期の心得 (広島…

3階のテンソルをL1ノルムの尺度で近似するランク１テンソルを見つけます，という論文を実装した．論文のタイトルにexact solutionがあるけど，L1ノルムの意味での最良１ランク近似をしてるわけではないので注意（かなり混乱した）． arxiv.org 最近，L1が流…

行列バランシングと言えば，sinkhornですが，今回はこの論文を実装しました． arxiv.org （論文では自然勾配法使ってますが，ここでは，１次までしか実装してないです） #include <iostream> #include <Eigen/Dense> #include <vector> #include <cmath> #include <time.h> #include <cstdio> #include <fstream> using names</fstream></cstdio></time.h></cmath></vector></eigen/dense></iostream>…

文字列hogehogeにhが何回出てくるか数えたい．ぐぐると str = "hogehoge" cnt = sum([1 for i = eachmatch(r"h", str)]) とかいいんじゃね？と書いてある． 強い人にはもっとシンプルに sum(split(str, "").=="h") で良いと教わった．なるほど． 追記 強い人…

3階のテンソルTをST-HOSVD, T-HOSVDより正確に高速に１ランク近似できる手法をJulia言語で実装した． あ，ここでいう階っていうのは，rank(T)じゃなくて，ndims(T)ね．3-waysとか，3-orderってこと．日本では，ndimsを階っていうので，（CP）ランクと混ざっ…

ついにDataFrameに入門してしまった...． 途中で列の名前を変えるにはこうする． using DataFrames names = ["colname1", "colname2", "colname3"] data = rand(3,3) df = DataFrame(data) rename!(df, Symbol.(names)) Symbolに変えないといけないことに気…

最近，毎日Juliaを書いているので，買って読んでみた．結論から言うと本当に良い買い物をしたと思う． 薄い技術書だけど、内容は評判通りかなり濃厚．Juliaの基本的な言語機能の他にも、コンパイラへのヒントの書き方、メモリ割り当てを削減するtips、プロフ…