12月第3週の振り返り - Napier

12月第3週（12/16-12/22）

債権総論

共同抵当とは、その語感とは異なり複数の債権者が抵当権を共有することではなく、同一の債権を担保するために複数の不動産に抵当権を設定することである。抵当権が実行された場合に後順位抵当権者の扱いを確定させるため、不動産間で負担の調整をする必要がある。

並列プログラミング

LU分解の過程において用いられるAの要素は求めようとするLまたはUの要素と同じ位置にある。そのためAの二次元配列とLUを組み合わせた二次元配列とを別々に持つ必要がないことが分かる（LまたはUの要素を求めたら、その位置にあるAの要素はもはや参照されることがないのでLまたはUで上書きしても支障がない）。これにより大規模な行列に対してメモリ使用量を削減することができる。

逐次計算によるLU分解は以下のようになる。素朴な実装ならLとUの配列をAとは別に確保しているところだが、省メモリ実装ではAを書き換えていくので、流れを追いやすいようプリプロセッサディレクティブを使ってL, Uと書けるようにした。

#include <stdio.h>
#include <math.h>
#define L A
#define U A
#define N 4

int main(void) {
    double A[N][N] = {
         0,  1,  4, -6,
         8,  1,  6, -2,
        -2,  0,  0,  3,
         4, -2, -2,  4
    };
    double b[N] = {-1, 13, 1, 4};
    
    double x[N];
    double c[N];

    // LU分解
    for (int k=0; k<N; k++) {
        // 部分ピボット選択
        // ピボットとする予定のU[k][k]が0だったり、小さかったりすると正しく計算できないので
        // 未処理部分の中で同列要素の絶対値が最大である行を探して入れ替える。対応するbも入れ替える
        // Aは途中までLとUに書き換えられているが、Lの未完成行は入れ替え可能なので行の全要素を入れ替えてよい
        double current_largest_abs = 0;
        int swap_row_idx;
        for(int i=k; i<N; i++){
            double candidate_abs = fabs(A[i][k]);
            if(current_largest_abs < candidate_abs){
                current_largest_abs = candidate_abs;
                swap_row_idx = i;
            }
        }
        // A swap
        for(int j=0; j<N; j++){
            double tmp = A[k][j];
            A[k][j] = A[swap_row_idx][j];
            A[swap_row_idx][j] = tmp;
        }
        // b swap
        double tmp = b[k];
        b[k] = b[swap_row_idx];
        b[swap_row_idx] = tmp;

        double pivot = A[k][k];
        double mult_inv = 1.0 / pivot; // 繰り返し除算をするのは遅いため逆数をとる

        // L 第k列(0-based index)
        for (int i=k+1; i<N; i++) {
            L[i][k] = A[i][k] * mult_inv;
        }


        // k+1行から最後の行までの全てにおいて
        // いつかは外側要素であるL[i][k]*U[k][j]を減算することになるので今のうちに行う
        for (int i=k+1; i<N; i++) {
            double Lik = L[i][k];
            
            // k+1列から最後の列までの全てにおいて減算する
            for (int j=k+1; j<N; j++) {
                U[i][j] = A[i][j] - Lik * U[k][j];
                // U 第k+1行(0-based index)はこれで求まる
                // それより下の行はまだ引くべき積が残っている
            }
        }
    }

    // Lc=bの上から順にcを求める
    for(int k=0; k<N; k++){
        c[k] = b[k];
        for(int j=0; j<k; j++){
            c[k] -= L[k][j] * c[j];
        }
    }

    // Ux=cの下から順にxを求める
    for(int k=N-1; k>=0; k--){
        x[k] = c[k];
        for(int j=N-1; j>k; j--){
            x[k] -= U[k][j] * x[j];
        }
        x[k] = x[k] / U[k][k];
    }

    // show x
    printf("x\n");
    for(int i=0; i<N; i++){
        printf("%.5lf ", x[i]);
    }
    return 0;
}

MPI並列化を施したものは授業の課題そのものになってしまうため載せない。

LU分解、前進消去、後退代入とそれぞれ手計算でやる分には迷うことは何もないのであるが、繰り返し同一の処理を行うプログラムとして表現されると直感から外れるので理解に手間取った。手計算ならば3x3行列で $l_{32}=\frac{a_{32}-l_{31}u_{12}}{u_{22}}$ と一気に要素の値を求めてしまうところ、プログラムでは $l_{32}$ に対する $-l_{31}u_{12}$ の処理を $u_{22}$ に対する $-l_{21}u_{12}$ と同一のタイミングで行っている。大きい行列の下の方になると、この減算が何回も積み重なり、最後にようやく除算（逆数を乗算）されて本当の値が求まる。

手計算では減算途中の値を書き留めるなどやっていられないので一発で要素を確定させる必要があり、減算回数が下に行くほど増えていくが、プログラムではメモリという完璧な記憶を活用して、最外ループが回るたび未確定の要素に対する段階的な減算を行っておくことができる。もちろん人間と同じやり方をしても計算の無駄は生じないのだが、PE間で行列を分散して保持している場合には $u_{ij}$ を求めるために複数のPEにまたがって $l_{ik},u_{kj}$ の値を送ってもらわねばならず、通信の無駄が生じる。

前進消去、後退代入においても同様に段階的な減算を行うが、行列ベクトル積の形なのでループごとの動きが容易に読めるだろう。

法医学

日本における身元不明死体の数は人口に比してかなり多い。また身元の確認ミスも起きている。顔貌による確認が主体であることに由来する。加えて高齢者の所在不明も多い。

乳歯・永久歯の萌出時期は決まっているため、生えそろい具合によって年齢を推定することができる。

個人識別の三大手段は指紋・歯科所見・DNA型であり、歯科所見は具体的に口腔内写真、レントゲン（デンタルX線）写真、デンタルチャートからなる。

CTによって口腔周辺の金属の発見ができ、インプラントや整形プレートなど身元特定の手がかりになることがある。

4歳程度までの歯は養育者のケアが必須であり、ケアがなければ齲歯（虫歯）になってしまう。これはデンタルネグレクトと呼ばれる。

ピンクの歯は自然の歯（生活歯）で、白っぽいのは根管治療を行って失活した歯である。X線撮影を行うと根管充填材は不透過なので明確に判定できる。

災害時には生前のカルテ・写真データを喪失してしまうことがありデータベース化が提言されている。

法遺伝学で行われるDNA型鑑定とはDNAの一部におけるバンドパターンなどを用いて特定に役立てるものである。Y染色体の多型（個体差）を利用して強姦事件の照合に使うこともあるが、Y染色体は親子・兄弟で共通しているため多型が稀ではなく、多くの男性で特徴が共通してしまうことがあり冤罪を生じさせうる。

科警研は以前はVNTR多型（MCT118）の鑑定であったが、現在はIdentifiler™を用いたSTR多型の鑑定に切り替わっている。同一の多型特徴を持つ二者が存在する確率は格段に下がったが、そもそも試料の採取・保存の過程におけるコンタミネーションなど取り扱いに起因する問題が生じる可能性は残っているため誤判の可能性を意識する必要がある。

現在では科捜研がDNA検査を事実上独占しており、捜査側による鑑定の公正性に懸念が生じている。

メディアプログラミング

tf-idfによる単語の重み付けを利用した特徴語抽出を行った。tf-idfは文書別に各語の出現頻度を集計し、全ての文書に広く出てくる単語はありふれたもの、少数の文書で集中的に現れる単語はそれらの文書を特徴づけるものと判断する。ただしこの手法は通常、新たに出会った文書がどの文書と類似しているかを判断することに用いられるので、出現頻度が少なすぎる語は判断に使えないとみなしてカットすることが多い。

例えばニュースサイトがジャンル分けなしに類似するニュースをリコメンドすることを考えたとき、「ネズミ」という語を特徴語として使えば、「ネズミ」という単語が含まれた新しいニュースに対しては動物や植物に関する過去のニュースをリコメンドすればよさそうだと判断することができるが、「アメリカミズハタネズミ」という特徴語は（そもそも名前の途中で分割されてしまいそうではあるが、分割されなかったとして）希少すぎるため新たに出会うニュースの中にまず出てこず利用チャンスがない。

また、たとえ頻度がある程度高くても1、2件のニュースの中でのみ何度も言及されているなら類似したニュースを見つけることには役立たないのでカットするということもありうる。

import numpy as np
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.metrics.pairwise import cosine_similarity

corpus = np.array(wakati_news_list)  # 単語分割済みのニュース群

vectorizer = TfidfVectorizer(max_features=1000, max_df=0.5, min_df=0.05)  # 5%以上半分以下のニュースで出現する語の頻度上位1000個についてtf-idfを計算する場合
corpus_tfidf = vectorizer.fit_transform(corpus)

target_tfidf = vectorizer.transform([wakati_new_news_text])  # 先に抽出した1000個の特徴語について新しいニュース文中での重要度をベクトルで出す

similarity_matrix = cosine_similarity(target_tfidf, corpus_tfidf)[0]  # ニュース群と新しいニュースのtf-idfベクトルの各コサイン類似度が得られる

授業とは関連しないが、以前ファイル出力でf.write()やprint(file=f)を使った際に一番良いファイル出力方法は何だろうかと思った。リストとジェネレータでメモリ使用量とコストを比較するによればジェネレータは巨大なデータ群を反復的に出力するとき格段にメモリ使用量が削減できることが示されている。これは自由に要素にアクセスできないことと引き換えに全ての要素をメモリに展開せずにおくからで、Effective Pythonでもジェネレータが推奨されている。

ただし多少遅くなる。速度を求める場合はメモリに載せられる限度で長いリストをstr.joinして一回でwrite()するのが良さそうだ。write()の複数呼び出しや文字列の逐次結合は遅くなる。日常的なキロバイト、メガバイトレベルの書き込みならばwrite(str.join(list))一発でやろうと思った。