[pukiwiki]
（誤差などを含んでおり、厳密には解くのが難しい）　連立一次方程式の近似解を、射影を使って解く方法が雑誌　「数学セミナー」 に掲載されていたので、　Pythonで実際に書いてみることにしました。
[/pukiwiki]

[pukiwiki]
参考文献
数学セミナー 2009年5月号
中村郁　「線形代数応用篇／たとえば歯の，CTスキャンのはなし 」
ERなどでもおなじみ、断層X線写真　CTスキャン。
人体の断面を輪切りの画像に変換するのに使われているアルゴリズムの解説記事。

具体的には、CTの画像を作成するのは、　（誤差などを含んでおり、厳密には解くのが難しい）　連立一次方程式の近似解を求めるのと同じ意味になるのだとか、、、　
-計測データを連立一次方程式だと考える
-一つ一つの一次方程式を　（超）平面として捉える。
-適当な点から　各面へ順番に射影を行っていくと、順番に　未知数 xへと収束していく

詳しくは　記事をご覧ください

数セミの記事の殆どは難しくて判らないのが多いのですが、この記事については、大変平易に書かれていて私でも読めました。

[[http://ecx.images-amazon.com/images/I/41t6ZffvmuL._SL160_.jpg:http://www.amazon.co.jp/exec/obidos/ASIN/B0021PPFPM/tamc-22/ref=nosim/]]

—-
以下コード

下記のような一次連立方程式の近似解を求めるプログラムです

a00*x0 + a01*x1 + a02*x2 + … +a0n*xn = b0
a10*x0 + a11*x1 + a12*x2 + … +a1n*xn = b1
…
ただし aは既知の係数 xは未知数

[/pukiwiki]

#psolve.py

from __future__ import division
from itertools import izip
from random import random
from math import sqrt
zip = izip
def iprod(a,b):
    "a and b must have same length"
    
    return sum(map( lambda (x,y):x*y,zip(a,b)))
def kprod(a,k):
    " a0*k,a1*k,,, an*k "
    return map(lambda x:x*k,a)
def vadd(a,b):
    return map( lambda (x,y):x+y, zip(a,b))
def vsub(a,b):
    return map( lambda (x,y):x-y, zip(a,b))
def nrm(a):
    return sqrt(sum([x**2 for x in a]))

def project(a,b,z):
    "project  z onto a0*x0+a1*x1+...an*xn=b"
    t=(b-iprod(a,z))/sum((x**2. for x in a))
    return vadd(z,kprod(a,t))
def psolve(alist,blist,z=None,loop=1):
    "pseudo solve"

    if loop>10 and z:
        print "@"
        return z
    #print "x",xlist
    
    if not z :
        z=[(random()-0.5)*100 for x in alist[0]]
    for i in xrange (500*loop):
        for a,b in zip(alist,blist):
            last=z
            z=project( a, b,z)
    if ( (nrm(z)/nrm(last))<0.9999):
        #print ".",last,z,".",
        
        z=psolve(alist,blist,z,loop=loop+1)
    return z

#ここからはテストコード
def iprodtest():
    def iprod2(a,b):
        ret=0
        for i,j in zip(a,b):
            ret+=i*j
        return ret
    a=[1,2,3]; b=[4,5,6]
    assert iprod(a,b)==iprod2(a,b)
    a=[100,110,1000]; b=[1000,0.1,0.5]
    assert iprod(a,b)==iprod2(a,b)
    from random import random
    for j in xrange(100):
        i=int(random()*100)
        a=[random() for x in xrange(i)]
        b=[random() for x in xrange(i)]
        assert iprod(a,b)==iprod2(a,b)

def psolve_test():
    xlist=[(random()-0.5)*100 for x in xrange(3)]

    alist=[[ (random()-0.5)*100 for x in xrange(3)] for y in xrange(3)]
    blist=[iprod(xlist,a) for a in alist]

    
    print u"与えられた問題　x",xlist
    z=psolve(alist,blist)
    a00,a01,a02=alist[0]
    a10,a11,a12=alist[1]
    a20,a21,a22=alist[2]
    print "行列式", ( -(a02*a11*a20) + a01*a12*a20 +
                   a02*a10*a21 - a00*a12*a21 -
                         a01*a10*a22 + a00*a11*a22 )
    print
    print u"未知数　z",z
    print u"誤差 e",vsub(xlist,z)
    #print alist
    assert nrm(vsub(xlist,z))< nrm(xlist)*0.2

if __name__=="__main__":
    iprodtest()
    psolve_test()
    print "done"

[pukiwiki]
上記テストコードだと、２０回に一回ぐらい誤差が大きいです。。。（解の収束が遅い？）
初期値や、テストデータの与え方が悪いのかしらん？
*一次方程式としての単純パーセプトロン
この関数の使い道として、こんなものを考えてみました。

一番シンプルなニューロコンピュータの仕組みに、「単純パーセプトロン」があります。
-[[ググる:単純パーセプトロン]]

単純パーセプトロンを数式として書くと

入力を　i0, i1,i2,,, in
それぞれの入力の重みを w0, w1, w2,,,wn　とした場合
w0*i0 + w1*i1 + w2*i2 + ... + wn*in = 出力値

となります。つまり、いくつかの入力値と出力値のサンプルがあれば、上記の、一次連立方程式の近似解を求める関数 psolveを使って、単純パーセプトロン　（シンプルなニューロコンピューター）のパラメータを求めることが出来る、、、はずです。

上記のテストコードでは近似解をそこそこ求めることが出来るようなので、また今度試してみたいと思います。
[/pukiwiki]