正方行列の固有値の実部が1未満なら行列の冪の極限は零行列

ジョルダン標準形を用いて示します。
ジョルダンブロック\(J(\lambda,k)\in\mathbb{R}^{k\times k}\)とは
{
J(\lambda,k)=\begin{pmatrix}
\lambda&1&0&\cdots&0&0\\
0&\lambda&1&0&\cdots&0\\
&&\ddots&&&\\
0&0&0&\cdots&\lambda&1\\
0&0&0&\cdots&0&\lambda
\end{pmatrix}
}
という形の行列です。
ジョルダンブロックの冪について\(n \geq k\)のときジョルダンブロックの冪は次の式で表せます
{
(J(\lambda,k))^n=\begin{pmatrix}
\lambda^n&\binom{n}{1}\lambda^{n-1}&\binom{n}{2}\lambda^{n-2}&\cdots&\binom{n}{k-2}\lambda^{n-k+2}&\binom{n}{k-1}\lambda^{n-k+1}\\
0&\lambda^n&\binom{n}{1}\lambda^{n-1}&\binom{n}{2}\lambda^{n-2}&\cdots&\binom{n}{k-2}\lambda^{n-k+2}\\
&&\ddots&&&\\
0&0&\cdots&\cdots&\lambda^n&\binom{n}{1}\lambda^{n-1}\\
0&0&0&\cdots&\cdots&\lambda^n
\end{pmatrix}
}
(証明)
帰納法と二項係数に関する関係式
{
\binom{n+1}{k}=\binom{n}{k}+\binom{n}{k-1}
}
から示すことができます。
全ての行列はジョルダンブロックの直和で表すことができるので(ジョルダン標準形を持つ)ので正方行列\(A\)の固有値の実部が全て1未満ならばジョルダンブロックの冪が収束して
{
\displaystyle\lim_{n\to\infty}A^n=O
}
が成り立ちます。

今日解いたぶんをはる。(100本ノック)

2章を解き始めた。

言語処理100本ノックをやってみる。

pythonに慣れていきたい。

www.cl.ecei.tohoku.ac.jp

こちらの教材が素晴らしいということだったので、初心者的に始めてみることにしました。

基本は文字列の操作についてなどなどが学べます。

最後の方はなかなか難しいらしいので、果たしてやり切れるのやら、、、。

以下は解いてみたものを貼り付けているだけです。

あと単語に区切るような操作は絶対正規表現使った方がいいんだろうなと思いながら、","とか":"とかの対応もせずにやってしまったので、
正規表現の章に入ったら、ちゃんとしていきたいですね。

#00. 文字列の逆順
#文字列"stressed"の文字を逆に(末尾から先頭に向かって)並べた文字列を得よ.

s="stressed"
#0   1   2   3   4   5   6   7   8
#| s | t | r | e | s | s | e | d |
#-9 -8  -7  -6  -5  -4  -3  -2  -1
#x[arg1:arg2:arg3]
#x[::-1]で逆順になる。
rev_s=s[::-1]
print(rev_s)
#00. 文字列の逆順
#文字列"stressed"の文字を逆に(末尾から先頭に向かって)並べた文字列を得よ.

s="stressed"
#01. 「パタトクカシーー」
#「パタトクカシーー」という文字列の1,3,5,7文字目を取り出して連結した文字列を得よ.
import io,sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
s=u"パタトクカシーー"
s1=""
for i in [1,3,5,7]:
    s1+=s[i-1]
print(s1)
#02. 「パトカー」+「タクシー」=「パタトクカシーー」
#「パトカー」+「タクシー」の文字を先頭から交互に連結して文字列「パタトクカシーー」を得よ.
import io,sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
s1=u"パトカー"
s2=u"タクシー"
s3=""
for (x,y) in zip(s1,s2):
    s3+=x+y
print(s3)
#03. 円周率
#"Now I need a drink, alcoholic of course, after the heavy lectures involving quantum mechanics."という文を
#単語に分解し,各単語の(アルファベットの)文字数を先頭から出現順に並べたリストを作成せよ.
import io,sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
s=u"Now I need a drink, alcoholic of course, after the heavy lectures involving quantum mechanics."
s=s.replace(",","")
s=s.replace(".","")
slist=s.split(" ")
l=[]
for i in slist:
    l.append(len(i))
print(l)
#04. 元素記号
#"Hi He Lied Because Boron Could Not Oxidize Fluorine. New Nations Might Also Sign Peace Security Clause. Arthur King Can."
#という文を単語に分解し,1, 5, 6, 7, 8, 9, 15, 16, 19番目の単語は先頭の1文字,それ以外の単語は先頭に2文字を取り出し,
#取り出した文字列から単語の位置(先頭から何番目の単語か)への連想配列(辞書型もしくはマップ型)を作成せよ.
import io,sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
s=u"Hi He Lied Because Boron Could Not Oxidize Fluorine. New Nations Might Also Sign Peace Security Clause. Arthur King Can."
s=s.replace(",","")
s=s.replace(".","")
slist=s.split(" ")
n=len(slist)
dic={}
one=[1,5,6,7,8,9,15,16,19]
for i in range(n):
    if i+1 in one:
        dic[slist[i][0:1]]=i+1
    else:
        dic[slist[i][0:2]]=i+1
print(dic)
#05. n-gram
#与えられたシーケンス(文字列やリストなど)からn-gramを作る関数を作成せよ.この関数を用い,"I am an NLPer"という文から単語bi-gram,文字bi-gramを得よ.
import io,sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
def word_ngram(s,n):
    s=s.replace(",","")
    s=s.replace(".","")
    sw=s.split(" ")
    ngram=[]
    for i in range(len(sw)-n+1):
        chank=""
        for j in range(i,i+n):
            if(j==i+n-1):
                chank+=sw[j]
            else:
                chank+=sw[j]+" "
        ngram.append(chank)
    return ngram

def char_ngram(s,n):
    s=s.replace(",","")
    s=s.replace(".","")
    sc=list(s.replace(" ",""))
    ngram=[]
    for i in range(len(sc)-n+1):
        chank=""
        for j in range(i,i+n):
            chank+=sc[j]
        ngram.append(chank)
    return ngram
#06. 集合
#"paraparaparadise"と"paragraph"に含まれる文字bi-gramの集合を,それぞれ, XとYとして求め,
#XとYの和集合,積集合,差集合を求めよ.さらに,'se'というbi-gramがXおよびYに含まれるかどうかを調べよ.
import io,sys
def word_ngram(s,n):
    s=s.replace(",","")
    s=s.replace(".","")
    sw=s.split(" ")
    ngram=[]
    for i in range(len(sw)-n+1):
        chank=""
        for j in range(i,i+n):
            if(j==i+n-1):
                chank+=sw[j]
            else:
                chank+=sw[j]+" "
        ngram.append(chank)
    return ngram

def char_ngram(s,n):
    s=s.replace(",","")
    s=s.replace(".","")
    sc=list(s.replace(" ",""))
    ngram=[]
    for i in range(len(sc)-n+1):
        chank=""
        for j in range(i,i+n):
            chank+=sc[j]
        ngram.append(chank)
    return ngram

sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
x="paraparaparadise"
y="paragraph"
X=set(char_ngram(x,2))
Y=set(char_ngram(y,2))
X_or_Y=X|Y
X_and_Y=X&Y
X_sub_Y=X-Y
if "se" in X:
    print("se is in X")
if "se" in Y:
    print("se is in Y")
#07. テンプレートによる文生成
#引数x, y, zを受け取り「x時のyはz」という文字列を返す関数を実装せよ.さらに,x=12, y="気温", z=22.4として,実行結果を確認せよ.
import io,sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
x=12
y="気温"
z="22.4"
print("%s時の%sは%s" % (x,y,z))
#08. 暗号文
#与えられた文字列の各文字を,以下の仕様で変換する関数cipherを実装せよ.

#英小文字ならば(219 - 文字コード)の文字に置換
#その他の文字はそのまま出力
#この関数を用い,英語のメッセージを暗号化・復号化せよ.
import io,sys
sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8')
sys.stderr = io.TextIOWrapper(sys.stderr.buffer, encoding='utf-8')
s="This is a 5000kg pen."

def encode(s):
    ss=list(s)
    sss=[]
    for i in ss:
        if i.islower():
            sss.append(chr(219-ord(i)))
        else:
            sss.append(i)
    s_code=""
    for i in sss:
        s_code+=i
    return s_code
def decode(s):
    ss=list(s)
    sss=[]
    for i in ss:
        if i.islower():
            sss.append(chr(219-ord(i)))
        else:
            sss.append(i)
    s_code=""
    for i in sss:
        s_code+=i
    return s_code

print(decode(encode(s)))
#09. Typoglycemia
#スペースで区切られた単語列に対して,各単語の先頭と末尾の文字は残し,
#それ以外の文字の順序をランダムに並び替えるプログラムを作成せよ.
#ただし,長さが4以下の単語は並び替えないこととする.適当な英語の文
#(例えば"I couldn't believe that I could actually understand what I was reading : the phenomenal power of the human mind .")
#を与え,その実行結果を確認せよ.
import random
def rand_change(s):
    sss=s.split(" ")
    sentence=[]
    for i in sss:
        if "," in i:
            ii.insert()
        if len(i)>4:
            ii=list(i[1:len(i)-1])
            random.shuffle(ii)
            ii.insert(0,i[0:1])
            ii.append(i[len(i)-1:len(i)])
            sentence.append("".join(ii))
        else:
            sentence.append(i)
    s_rand=" ".join(sentence)
    return s_rand
s="I couldn't believe that I could actually understand what I was reading : the phenomenal power of the human mind ."
print(rand_change(s))

ネットビジネス系のアカウントを半自動ブロックするやつ書いた。

タイトルの通りです。
ネットビジネスとか、アフィリエイトとか流行ってますよね。
非常にくさいです。
自動でブロックできないかなーって思って探してたら、pythonでめっちゃお手軽にかけるライブラリを発見。

ここの記事をそのまま参考にさせていただいきました。

kivantium.hateblo.jp


今は完全に自分用だけども気が向いたら、アプリ的などにできるといいですね。

コードも乗っけておきます。アクセストークンのとことnameのとこを自分のものに変えれば、誰でも使えるとは思います。

#-*- coding: utf-8 -*-
import tweepy
#自分のアクセストークンを持って来る
import settings_m as settings
#認証を行う
auth = tweepy.OAuthHandler(settings.CK, settings.CS)
auth.set_access_token(settings.AT, settings.AS)
api = tweepy.API(auth)
#指定したユーザーのフォロワーをリストで取得
followers=api.followers_ids("MAEA_2")
#ブロックするキーワードを指定
block_words=["ネットビジネス","稼げる","収益","アフィリエイト","マーケティング","好きな事をしながら","スマホ一つで","スマホ1つで"]
block=-1
#探索
for follower in followers:
    block=0
    user=api.get_user(follower)
    for word in block_words:
        #ユーザー名かプロフィールにブロックワードが含まれていたらブロック
        if word in user.name or word in user.description:
            block=1
    if block==1:
        print(user.screen_name)
        api.create_block(follower)
    print(block)

追記:アクセス制限があるので、フォロワーが多いと最後まで回らないかもしれない。

a^2/(b+c)+b^2/(c+a)+c^2/(a+b)が自然数となる自然数の組について 衝撃の解答編

この前の記事で\(a,b,c\in\mathbb{N}\)で
\( a < b < c \)とするとき
\[\frac{a^2}{b+c}+\frac{b^2}{c+a}+\frac{c^2}{a+b}\in\mathbb{N}\]
となるような組が沢山ある(おそらく無限?)ってことは前の記事で書いたんですが
さらに制限をつけて
\(a,b,c\)が互いに素
としたら解として
\[a=13,b=77,c=183\]
しかでてこなくて困ってしまい、math over flow に投げました。
数学の知恵袋サイト的なところで、日夜多くの数学の議論で賑わっています。

そうすると数時間後にはすぐに回答をいただくことができました!!

mathoverflow.net

回答者様の投稿によると\(a+b+c\)の順に解を並べたとき、次に小さな解で互いに素なものは


\(a=15349474555424019\)
\(b=35633837601183731\)
\(c=105699057106239769\)


だそうです。これは流石に全探索では見つけられない...

ちなみにこのとき
\[\frac{a^2}{b+c}+\frac{b^2}{c+a}+\frac{c^2}{a+b}=231293021292774912\]

となります。ひょえええってかんじですね。

背後に理論があるようでそのことについても書いてくださっており、非常に勉強になりました!

追記:
さらに小さな解についても見つけてくださってますね。
\[a=248,b=2755,c=7227\]


\[\frac{a^2}{b+c}+\frac{b^2}{c+a}+\frac{c^2}{a+b}=18414\]
となり、確かに解です!

a>b>cでa^2/(b+c)+b^2/(c+a)+c^2/(a+b)が自然数となる組について

\(a > b>c\)となる自然数の組のうち
\[\frac{a^2}{(b+c)}+\frac{b^2}{(c+a)}+\frac{c^2}{(a+b)}\]
自然数となるようなものを考えたい。

現在見つかっている組
\((a,b,c)=\)
\((3k,2k,1k) \min k=5\)
\((5k,2k,1k) \min k=7\)
\((7k,3k,2k) \min k=5\)
\((5k,4k,3k) \min k=7\)
\((4k,2k,1k) \min k=10\)
\((5k,3k,1k) \min k=8\)
\((5k,4k,1k) \min k=9\)
\((11k,4k,1k) \min k=5\)
\((7k,5k,1k) \min k=8\)
\((7k,5k,3k) \min k=8\)
\((7k,5k,2k) \min k=9\)
\((33k,7k,2k) \min k=2\)
\((23k,7k,5k) \min k=3\)
\((5k,3k,2k) \min k=14\)
\((17k,3k,1k) \min k=5\)
\((17k,11k,4k) \min k=5\)
\((8k,3k,1k) \min k=11\)
\((8k,7k,3k) \min k=11\)
\((11k,5k,4k) \min k=9\)
\((3k,2k,1k) \min k=5\)
\((5k,2k,1k) \min k=7\)
\((7k,3k,2k) \min k=5\)
\((5k,4k,3k) \min k=7\)
\((4k,2k,1k) \min k=10\)
\((5k,3k,1k) \min k=8\)
\((5k,4k,1k) \min k=9\)
\((11k,4k,1k) \min k=5\)
\((7k,5k,1k) \min k=8\)
\((7k,5k,3k) \min k=8\)
\((7k,5k,2k) \min k=9\)
\((33k,7k,2k) \min k=2\)
\((23k,7k,5k) \min k=3\)
\((5k,3k,2k) \min k=14\)
\((17k,3k,1k) \min k=5\)
\((17k,11k,4k) \min k=5\)
\((8k,3k,1k) \min k=11\)
\((8k,7k,3k) \min k=11\)
\((11k,5k,4k) \min k=9\)
\((8k,7k,5k) \min k=13\)
\((11k,7k,3k) \min k=10\)
\((11k,9k,1k) \min k=10\)
\((16k,9k,5k) \min k=7\)
\((8k,6k,1k) \min k=14\)
\((14k,11k,10k) \min k=8\)
\((19k,11k,1k) \min k=6\)
\((23k,7k,1k) \min k=5\)
\((7k,2k,1k) \min k=18\)
\((9k,5k,1k) \min k=14\)
\((22k,5k,3k) \min k=6\)
\((6k,5k,4k) \min k=22\)
\((19k,16k,5k) \min k=7\)
\((17k,13k,7k) \min k=8\)
\((4k,3k,1k) \min k=35\)
\((9k,7k,5k) \min k=16\)
\((8k,7k,2k) \min k=18\)
\((29k,27k,13k) \min k=5\)
\((17k,7k,1k) \min k=9\)
\((11k,3k,1k) \min k=14\)
\((23k,5k,2k) \min k=7\)
\((33k,23k,7k) \min k=5\)
\((11k,5k,1k) \min k=16\)
\((13k,8k,2k) \min k=14\)
\((7k,6k,1k) \min k=26\)
\((183k,77k,13k) \min k=1\)
\((93k,7k,5k) \min k=2\)
\((17k,7k,4k) \min k=11\)
\((27k,13k,8k) \min k=7\)
\((9k,6k,1k) \min k=21\)
\((13k,8k,7k) \min k=15\)
\((41k,8k,7k) \min k=5\)
\((23k,13k,7k) \min k=9\)
\((13k,5k,3k) \min k=16\)
\((19k,5k,3k) \min k=11\)
\((19k,9k,2k) \min k=11\)
\((13k,12k,5k) \min k=17\)
\((23k,17k,7k) \min k=10\)
\((11k,10k,1k) \min k=21\)
\((13k,7k,5k) \min k=18\)
\((49k,17k,11k) \min k=5\)
\((129k,25k,11k) \min k=2\)
\((13k,3k,2k) \min k=20\)
\((65k,33k,7k) \min k=4\)
\((53k,13k,2k) \min k=5\)
\((19k,5k,2k) \min k=14\)
\((19k,11k,9k) \min k=14\)
\((9k,3k,1k) \min k=30\)
\((34k,21k,11k) \min k=8\)
\((17k,13k,3k) \min k=16\)
\((34k,29k,6k) \min k=8\)
\((21k,5k,4k) \min k=13\)
\((14k,6k,1k) \min k=20\)
\((7k,3k,1k) \min k=40\)
\((4k,3k,2k) \min k=70\)
\((15k,4k,1k) \min k=19\)
\((43k,37k,5k) \min k=7\)
\((19k,11k,5k) \min k=16\)
\((28k,27k,17k) \min k=11\)
\((14k,9k,1k) \min k=23\)
\((17k,3k,2k) \min k=19\)
\((25k,14k,1k) \min k=13\)
\((15k,13k,7k) \min k=22\)
\((6k,2k,1k) \min k=56\)
\((15k,8k,1k) \min k=23\)
\((15k,13k,8k) \min k=23\)
\((25k,11k,3k) \min k=14\)
\((25k,17k,3k) \min k=14\)
\((14k,11k,1k) \min k=25\)
\((71k,13k,7k) \min k=5\)
\((51k,29k,19k) \min k=7\)
\((8k,7k,1k) \min k=45\)
\((28k,11k,5k) \min k=13\)
\((23k,13k,3k) \min k=16\)
\((11k,6k,4k) \min k=34\)
\((6k,3k,1k) \min k=63\)
\((11k,4k,3k) \min k=35\)
\((35k,19k,9k) \min k=11\)
\((7k,4k,1k) \min k=55\)
\((7k,4k,3k) \min k=55\)
\((15k,11k,7k) \min k=26\)
\((13k,11k,2k) \min k=30\)
\((17k,15k,8k) \min k=23\)
\((9k,7k,2k) \min k=44\)
\((16k,11k,9k) \min k=25\)
\((31k,21k,5k) \min k=13\)
\((37k,7k,5k) \min k=11\)
\((18k,7k,5k) \min k=23\)
\((32k,17k,7k) \min k=13\)
\((19k,6k,5k) \min k=22\)
\((6k,4k,1k) \min k=70\)
\((25k,11k,9k) \min k=17\)
\((11k,2k,1k) \min k=39\)
\((11k,7k,2k) \min k=39\)
\((29k,11k,4k) \min k=15\)
\((19k,4k,1k) \min k=23\)
\((20k,19k,13k) \min k=22\)
\((37k,19k,5k) \min k=12\)
\((25k,19k,11k) \min k=18\)
\((41k,25k,19k) \min k=11\)
\((13k,2k,1k) \min k=35\)
\((35k,33k,17k) \min k=13\)
\((19k,5k,1k) \min k=24\)
\((12k,7k,2k) \min k=38\)
\((27k,8k,7k) \min k=17\)
\((17k,11k,10k) \min k=27\)
\((6k,5k,1k) \min k=77\)
\((11k,10k,4k) \min k=42\)
\((235k,17k,3k) \min k=2\)
\((19k,14k,11k) \min k=25\)
\((14k,4k,3k) \min k=34\)
\((9k,2k,1k) \min k=55\)
\((11k,7k,4k) \min k=45\)
\((31k,25k,17k) \min k=16\)
\((18k,17k,10k) \min k=28\)
\((46k,35k,9k) \min k=11\)
\((37k,15k,13k) \min k=14\)
\((10k,3k,2k) \min k=52\)
\((13k,7k,1k) \min k=40\)
\((13k,11k,9k) \min k=40\)
\((20k,19k,6k) \min k=26\)
\((25k,3k,2k) \min k=21\)
\((59k,31k,4k) \min k=9\)
\((77k,63k,58k) \min k=7\)
\((20k,7k,1k) \min k=27\)
\((17k,15k,1k) \min k=32\)
\((21k,19k,5k) \min k=26\)
\((43k,22k,17k) \min k=13\)
\((44k,21k,5k) \min k=13\)
\((41k,29k,1k) \min k=14\)
\((20k,9k,7k) \min k=29\)
\((83k,13k,8k) \min k=7\)
\((53k,37k,35k) \min k=11\)
\((9k,4k,1k) \min k=65\)
\((13k,7k,2k) \min k=45\)
\((31k,14k,5k) \min k=19\)
\((18k,15k,7k) \min k=33\)
\((17k,4k,3k) \min k=35\)
\((35k,33k,1k) \min k=17\)
\((35k,33k,16k) \min k=17\)
\((23k,9k,4k) \min k=26\)
\((86k,19k,9k) \min k=7\)
\((43k,22k,13k) \min k=14\)
\((55k,33k,17k) \min k=11\)
\((19k,17k,15k) \min k=32\)
\((12k,5k,3k) \min k=51\)
\((28k,5k,2k) \min k=22\)
\((77k,67k,3k) \min k=8\)
\((28k,27k,5k) \min k=22\)
\((39k,11k,5k) \min k=16\)
\((57k,20k,13k) \min k=11\)
\((19k,8k,3k) \min k=33\)
\((5k,4k,2k) \min k=126\)
\((29k,15k,7k) \min k=22\)
\((31k,11k,4k) \min k=21\)
\((41k,7k,1k) \min k=16\)
\((41k,15k,1k) \min k=16\)
\((17k,10k,3k) \min k=39\)
\((35k,3k,1k) \min k=19\)
\((35k,22k,3k) \min k=19\)
\((24k,11k,4k) \min k=28\)
\((21k,11k,1k) \min k=32\)
\((17k,7k,3k) \min k=40\)
\((20k,14k,1k) \min k=34\)
\((12k,7k,3k) \min k=57\)
\((38k,28k,17k) \min k=18\)
\((12k,9k,7k) \min k=57\)
\((43k,35k,13k) \min k=16\)
\((9k,5k,2k) \min k=77\)
\((7k,5k,4k) \min k=99\)
\((63k,47k,25k) \min k=11\)
\((24k,5k,1k) \min k=29\)
\((20k,15k,1k) \min k=35\)
\((9k,6k,4k) \min k=78\)
\((177k,112k,68k) \min k=4\)
\((21k,13k,5k) \min k=34\)
\((13k,9k,2k) \min k=55\)
\((11k,9k,2k) \min k=65\)
\((11k,9k,4k) \min k=65\)
\((55k,49k,1k) \min k=13\)
\((8k,2k,1k) \min k=90\)
\((20k,8k,1k) \min k=36\)
\((15k,9k,1k) \min k=48\)
\((9k,7k,1k) \min k=80\)
\((103k,23k,12k) \min k=7\)
\((7k,6k,2k) \min k=104\)
\((35k,21k,19k) \min k=21\)
\((44k,7k,5k) \min k=17\)
\((11k,6k,5k) \min k=68\)
\((47k,9k,7k) \min k=16\)
\((20k,18k,7k) \min k=38\)
\((51k,19k,9k) \min k=15\)
\((9k,8k,1k) \min k=85\)
\((9k,8k,7k) \min k=85\)
\((35k,9k,1k) \min k=22\)
\((43k,29k,13k) \min k=18\)
\((25k,24k,7k) \min k=31\)
\((37k,23k,5k) \min k=21\)
\((21k,19k,16k) \min k=37\)
\((46k,29k,5k) \min k=17\)
\((29k,20k,7k) \min k=27\)
\((11k,7k,1k) \min k=72\)
\((159k,28k,17k) \min k=5\)
\((14k,5k,1k) \min k=57\)
\((25k,7k,3k) \min k=32\)
\((23k,22k,13k) \min k=35\)
\((26k,9k,5k) \min k=31\)
\((101k,31k,9k) \min k=8\)
\((21k,5k,3k) \min k=39\)
\((13k,5k,2k) \min k=63\)
\((9k,5k,4k) \min k=91\)
\((13k,8k,1k) \min k=63\)
\((13k,8k,5k) \min k=63\)
\((25k,8k,7k) \min k=33\)
\((49k,19k,2k) \min k=17\)
\((15k,13k,11k) \min k=56\)
\((21k,19k,17k) \min k=40\)
\((47k,13k,5k) \min k=18\)
\((10k,7k,5k) \min k=85\)
\((15k,13k,6k) \min k=57\)
\((19k,17k,1k) \min k=45\)
\((27k,13k,5k) \min k=32\)
\((67k,11k,10k) \min k=13\)
\((46k,11k,9k) \min k=19\)
\((8k,3k,2k) \min k=110\)
\((11k,9k,7k) \min k=80\)
\((23k,16k,5k) \min k=39\)
\((18k,17k,7k) \min k=50\)
\((26k,9k,1k) \min k=35\)
\((26k,19k,9k) \min k=35\)
\((38k,25k,7k) \min k=24\)
\((17k,10k,7k) \min k=54\)
\((17k,10k,8k) \min k=54\)
\((22k,20k,13k) \min k=42\)
\((37k,13k,2k) \min k=25\)
\((24k,15k,1k) \min k=39\)
\((8k,5k,1k) \min k=117\)
\((27k,8k,1k) \min k=35\)
\((35k,19k,1k) \min k=27\)
\((73k,57k,47k) \min k=13\)
\((481k,329k,9k) \min k=2\)
\((23k,13k,5k) \min k=42\)
\((13k,12k,3k) \min k=75\)
\((14k,7k,1k) \min k=70\)
\((29k,5k,1k) \min k=34\)
\((76k,41k,15k) \min k=13\)
\((247k,185k,153k) \min k=4\)
\((43k,3k,2k) \min k=23\)
\((31k,9k,1k) \min k=32\)

スパースモデリングを超初歩から始めてみた。(ISTAを実装した)

ブログを書くのかなり久しぶりな感じがしますね。

春休み、ちょっと勉強しようと思って東北大の大関正之先生の書かれたスパースモデリングの入門テキスト

「今日からできるスパースモデリング」 大関正之

を読んで、すごいおもしろいなーと感じたので、ちょっとやってみることにしました。

次のような問題を考えます。\(A\)を\(M\times N\)行列として\(N\)次元ベクトル\(\boldsymbol{x}\)についての線形連立方程式
\[\boldsymbol{y}=A\boldsymbol{x}\]
を考えます。

このとき\(N < M\)であれば与えられた方程式の数の方が多いので方程式は解けます。これを優決定系といいます。

これに対して\(N>M\)であるときは条件式が足りないので方程式は解けないということになりますが、

ここで\(\boldsymbol{x}\)に次のような条件が事前に与えられているとしましょう。

\(K\)スパース性:\(\boldsymbol{x}\)のほとんどの成分は\(0\)であり,非ゼロ成分の個数は\(K\)個である。

\(K< M\)であれば方程式を解くことができます。

しかしながら現実の場面においては非ゼロ成分がどこに分布してるのかということや\(K\)自体がわからない場合が多いです。

このような場合、次のような最適化問題を解くことになりますが、これはNP困難な問題となり計算量爆発を起こしてしまいます。

\[\min_{\boldsymbol{x}}\|x\|_0\,\,\, \mathrm{s.t.}\,\, \boldsymbol{y}=A\boldsymbol{x}\]

ここで、\(\|\cdot\|_0\)は\(L_0\)ノルムと呼ばれ、非ゼロ成分の個数を表します。

そのような場合に有効なアプローチとして圧縮センシングの方法があるそうです。

これは\(L_0\)ノルムの代わりに\(L_1\)ノルム\(\|\cdot\|_1\)を用いるものです。\(L_1\)ノルムとは

\[\|\boldsymbol{x}\|=|x_1|+|x_2|+\ldots+|x_N|\]

で定義されるノルムのことです。

\(L_1\)ノルムを用いて次の最適化問題を考えます。

\(L_1\)再構成: \(\min_{\boldsymbol{x}}\|x\|_1\,\,\, \mathrm{s.t.}\,\, \boldsymbol{y}=A\boldsymbol{x}\)


\(L_1\)ノルムは非ゼロ成分の個数が少ないほど小さくなりますが、成分の大きさを反映しているため\(L_0\)ノルムによる最適化とは本質的には異なるはずですが、多くの場合うまくいくらしいです。

\(L_1\)再構成を具体的に解く計算アルゴリズムの基本的なものにISTA(Iterative Shrinkage Soft-thresholding Algorithm)というのがあります。次のようなアルゴリズムです。
ISTA

  1. \(t=0\),\(\boldsymbol{x} [0]\)を適当に初期化する。例えば\(\boldsymbol{x} [0]=A^\top \boldsymbol{y}\)。
  2. \(\boldsymbol{v}[t]=x[t]+(L\lambda)^{-1}A^\top(\boldsymbol{y}-A\boldsymbol{x}[t])\)
  3. \(\boldsymbol{x}[t+1]=S_{L^{-1}}(\boldsymbol{v}[t])\),\(t=t+1\)
  4. 終了条件を満たすまで2.ー3.を繰り返す。

ISTAの導出などについては上記の資料を見てください。\(S_\lambda(\cdot)\)というのは軟判定閾値関数とよばれるもので

\[
S_\lambda(y)=
\begin{cases}
y-\lambda &&y>\lambda\\
0 &&-\lambda\leq y\leq\lambda\\
y+\lambda &&y<-\lambda
\end{cases}
\]

で定義されます。\(S_\lambda(\boldsymbol{y})\)は\(\boldsymbol{y}\)の各成分に\(S_\lambda\)を作用させたものという意味です。

以下がコードです。なんだかあまり上手く原情報ベクトルが再現されてなかったのでどこかで間違ってそうです。

間違いなど、ドシドシ教えてください。それではまた。

import numpy as np
import random

#軟判定閾値関数
def ST(y,lam):
    if y>lam :
        return y-lam
    elif y<-lam :
        return y+lam
    else:
        return 0
N=1000
M=100
lam=1E-10
ep=1E-8
#原情報ベクトルの非零成分の個数
K=20
# 観測行列を生成
A=np.random.randn(M,N)
A_T=A.transpose()
#リプシッッツ定数の設定 L=|A^TA|_F(フロベニウスノベル)
L=np.linalg.norm(A_T.dot(A),ord="fro")/lam
print("L=%f" % L)
# スパースな原信号ベクトルを生成
x_0=np.zeros(N)
r_K=random.sample(range(N),K)
for i in r_K:
    x_0[i]=np.random.rand()
#観測ベクトルを生成
y=A.dot(x_0)
#ISTAの初期値の設定
t=0
x_t=A_T.dot(y)
#収束するまで繰り返す
while np.linalg.norm(y-A.dot(x_t)) > ep:
    v_t=x_t+A_T.dot(y-A.dot(x_t))/L/lam
    print("|y-Ax[%d]|_2=%f" % (t,np.linalg.norm(y-A.dot(x_t))))
    for i in range(N):
        x_t[i]=ST(v_t[i],1/L)
    t=t+1
print(np.linalg.norm(x_t-x_0))