qr с выбором

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
__author__ = 'Lenovo'
import pprint
from pprint import pprint
from math import sqrt
import numpy as np
np.set_printoptions(suppress=True, linewidth=160, precision=10)
import scipy
from scipy import linalg
import pprint
A = scipy.array([[1,2,3], [1,5,6], [1,8,9], [1,11,12]])
Q, R, P = scipy.linalg.qr(A, pivoting=True)
print "A:"
pprint.pprint(A)
print "Q:"
pprint.pprint(Q)
print "R:"
pprint.pprint(R)
print "P:"
pprint.pprint(P)
####################################################################################
def sign(x):
return cmp(x, 0)
def norm(x):
return sum(i ** 2 for i in x) ** 0.5
def house(x):
n = len(x)
mu = norm(x)
v = x
if mu != 0:
beta = x[0] + sign(x[0]) * mu
v = v / beta
v[0] = 1
return v
def rowHouse(a, v):
m = len(v)
beta = -2.0 / sum(i ** 2 for i in v)
vm = np.array([v]).transpose()
w = beta * np.dot(a.transpose(), vm)
p = np.dot(vm, w.transpose())
pa = a + p
return pa
def swap(a, r, k):
temp = a[:, k].copy()
a[:, k] = a[:, r]
a[:, r] = temp
def swapElems(c, r, k):
temp = c[k]
c[k] = c[r]
c[r] = temp
def qrColumnPivoting(permutation, a):
m = len(a)
n = len(a[0])
c = np.zeros(n)
for i in range(0, n):
c[i] = sum(j**2 for j in a[:,i])
r = 0
t = c.max()
k = c.argmax()
while t > 0:
permutation[r] = k
permutation[k] = r
swap(a, r, k)
swapElems(c, r, k)
v = house(a[r:m, r])
a[r:m, r:n] = rowHouse(a[r:m, r:n], v)
a[r+1:m, r] = v[1:m]
for i in range(r+1, n):
c[i] -= a[r, i]**2
if r < n-1:
t = c[r+1:].max()
k = c[r+1:].argmax()
else:
t = 0
r += 1
return permutation, a, r
def fod(a):
return qt, u
a = np.array([[1,2,3], [1,5,6], [1,8,9], [1,11,12]], dtype=float)
p, a, rank = qrColumnPivoting(np.arange(len(a[0])), a)
qt, u = fod(a)
print "A:"
pprint.pprint(a)
print "P:"
pprint.pprint(p)