Sayısal olarak ortogonal projeksiyon

Question

Numpy.linalg.lstsq - yöntemiyle bir düzlem hesapladığım 3B noktaların bir listesi var. Ama Şimdi bu düzlem içine her bir nokta için bir ortogonal projeksiyon yapmak istiyorum, ama benim hata bulamıyorum:

from numpy.linalg import lstsq 

def VecProduct(vek1, vek2): 
    return (vek1[0]*vek2[0] + vek1[1]*vek2[1] + vek1[2]*vek2[2]) 

def CalcPlane(x, y, z): 
    # x, y and z are given in lists 
    n = len(x) 
    sum_x = sum_y = sum_z = sum_xx = sum_yy = sum_xy = sum_xz = sum_yz = 0 
    for i in range(n): 
     sum_x += x[i] 
     sum_y += y[i] 
     sum_z += z[i] 
     sum_xx += x[i]*x[i] 
     sum_yy += y[i]*y[i] 
     sum_xy += x[i]*y[i] 
     sum_xz += x[i]*z[i] 
     sum_yz += y[i]*z[i] 

    M = ([sum_xx, sum_xy, sum_x], [sum_xy, sum_yy, sum_y], [sum_x, sum_y, n]) 
    b = (sum_xz, sum_yz, sum_z) 

    a,b,c = lstsq(M, b)[0] 

    ''' 
    z = a*x + b*y + c 
    a*x = z - b*y - c 
    x = -(b/a)*y + (1/a)*z - c/a 
    ''' 

    r0 = [-c/a, 
      0, 
      0] 

    u = [-b/a, 
     1, 
     0] 

    v = [1/a, 
     0, 
     1] 

    xn = [] 
    yn = [] 
    zn = [] 

    # orthogonalize u and v with Gram-Schmidt to get u and w 

    uu = VecProduct(u, u) 
    vu = VecProduct(v, u) 
    fak0 = vu/uu 
    erg0 = [val*fak0 for val in u] 
    w = [v[0]-erg0[0], 
     v[1]-erg0[1], 
     v[2]-erg0[2]] 
    ww = VecProduct(w, w) 

    # P_new = ((x*u)/(u*u))*u + ((x*w)/(w*w))*w 
    for i in range(len(x)): 
     xu = VecProduct([x[i], y[i], z[i]], u) 
     xw = VecProduct([x[i], y[i], z[i]], w) 
     fak1 = xu/uu 
     fak2 = xw/ww 
     erg1 = [val*fak1 for val in u] 
     erg2 = [val*fak2 for val in w] 
     erg = [erg1[0]+erg2[0], erg1[1]+erg2[1], erg1[2]+erg2[2]] 
     erg[0] += r0[0] 
     xn.append(erg[0]) 
     yn.append(erg[1]) 
     zn.append(erg[2]) 

    return (xn,yn,zn) 

Bu bana bir düzlemde hepsi noktaların bir listesini döndürür, ama ben görüntülediğinizde Onlar, olması gereken pozisyonlarda değiller. Ben zaten bu sorunu çözmek için belirli bir yerleşik bir yöntem olduğuna inanıyorum, ama bulamadım herhangi bir = (

Answer 1

Çok düşük bir np.lstsq kullanımı yapıyorsunuz, çünkü önceden yapılandırılmış bir 3x3 matrisini besliyorsunuz, . bunun yerine işi yapmak icar böyle yapacağını:

import numpy as np 

def calc_plane(x, y, z): 
    a = np.column_stack((x, y, np.ones_like(x))) 
    return np.linalg.lstsq(a, z)[0] 

>>> x = np.random.rand(1000) 
>>> y = np.random.rand(1000) 
>>> z = 4*x + 5*y + 7 + np.random.rand(1000)*.1 
>>> calc_plane(x, y, z) 
array([ 3.99795126, 5.00233364, 7.05007326]) 

o z, yani kullanım sıfır olmama katsayısı bağlı değildir uçaktan için bir formül kullanmak aslında daha uygundur a*x + b*y + c*z = 1 Benzer şekilde a, b ve c hesaplamalarını da yapabilirsiniz:

Noktaları bir düzleme yansıtmak için, alternatif denklemi kullanarak, (a, b, c) vektörü düzleme diktir. (a, b, c)/(a**2+b**2+c**2) no.lu noktanın düzlemde olduğunu kontrol etmek kolaydır, bu nedenle tüm noktaları nokta üzerinde normal vektöre yansıtarak, o projeksiyonu noktalardan çıkararak, sonra da bu noktaya geri göndererek, projeksiyon, tüm noktaları düzlemde referans alarak yapılabilir. köken. Yani şimdi

def project_points(x, y, z, a, b, c): 
    """ 
    Projects the points with coordinates x, y, z onto the plane 
    defined by a*x + b*y + c*z = 1 
    """ 
    vector_norm = a*a + b*b + c*c 
    normal_vector = np.array([a, b, c])/np.sqrt(vector_norm) 
    point_in_plane = np.array([a, b, c])/vector_norm 

    points = np.column_stack((x, y, z)) 
    points_from_point_in_plane = points - point_in_plane 
    proj_onto_normal_vector = np.dot(points_from_point_in_plane, 
            normal_vector) 
    proj_onto_plane = (points_from_point_in_plane - 
         proj_onto_normal_vector[:, None]*normal_vector) 

    return point_in_plane + proj_onto_plane 

gibi bir şey yapabilirsiniz: Sen olarak aşağıdaki bunu yapabilir Sadece matrisleri her şeyin bir seçenektir yapabilirsiniz

>>> project_points(x, y, z, *calc_plane_bis(x, y, z)) 
array([[ 0.13138012, 0.76009389, 11.37555123], 
     [ 0.71096929, 0.68711773, 13.32843506], 
     [ 0.14889398, 0.74404116, 11.36534936], 
     ..., 
     [ 0.85975642, 0.4827624 , 12.90197969], 
     [ 0.48364383, 0.2963717 , 10.46636903], 
     [ 0.81596472, 0.45273681, 12.57679188]]) 

Answer 2

.

Eğer bir matris X için satır vektörleri olarak puanlarınızı eklerseniz ve y sonra bir vektör, en küçük kareler çözümü için parametreler vektörü beta olduğunu şunlardır:

import numpy as np 

beta = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y)) 

ama daha kolay bir yolu varsa, orada biz Projeksiyon yapmak istediğiniz: QR ayrışımı bize Q.T olarak bir ortonormal projeksiyon matrisi verir ve Q'un kendisi de ortonormal temel vektörlerinin matrisidir. Bu yüzden, önce QR'u oluşturabilir, ardından beta'u alabilir, ardından noktaları yansıtmak için Q.T'u kullanabiliriz.

QR:

Q, R = np.linalg.qr(X) 

beta:

# use R to solve for beta 
# R is upper triangular, so can use triangular solver: 
beta = scipy.solve_triangular(R, Q.T.dot(y)) 

Yani şimdi beta var ve çok basit Q.T kullanarak noktaları yansıtabilirsiniz:

X_proj = Q.T.dot(X) 

sne o!

Daha fazla bilgi ve grafiksel piccies falan istiyorsan

, ben de benzer bir şey yapıyor iken, notlar bir sürü yaptı: https://github.com/hughperkins/selfstudy-IBP/blob/9dedfbb93f4320ac1bfef60db089ae0dba5e79f6/test_bases.ipynb

(Düzenleme:, bir önyargı terimini eklemek istiyorsanız unutmayın böylece En iyi uyum kaynağından geçmek zorunda değilsiniz, sadece 1'leri içeren ek bir sütun ekleyebilmeniz için X ile bias terimi/özelliği gibi davranır)