This Gist demonstrates some code for generating Householder transformations, and using it to perform QR decompositions, including pytest unit tests. This is prototype code, and could be made more efficient in various ways; for example, when performing the QR decomposition, calculating v (the Householder vector) could be made a lot more efficient by using the fact that y is always a one-hot vector, and the memory efficiency can be improved by normalising v such that the first non-zero element is unity, and storing the rest of the vector in the corresponding column of R, as described in chapter 5 of Matrix Computations by Golub and Van Loan (4th edition).
import numpy as np
class Householder:
def __init__(self, x, y):
"""
Initialise a householder transformation which transforms the vector x to
point in the direction of a vector y, by reflecting x in a hyperplane
through the origin, whose normal is denoted by v. The matrix of this
linear transformation is denoted by P and is equal to I -
2*[email protected]/(v.T@v) == I - beta*[email protected], where beta = 2/(v.T@v). By
definition, P@x == alpha*y for some scalar alpha.
According to the derivation in section 5.1.2 of Matrix Computations by
Golub and Van Loan (4th edition), expanding P = I - 2*[email protected]/(v.T@v) and
P@x == alpha*y implies that v is parallel to x - alpha*y, and the
magnitude of v is not unique, so we can set v = x - alpha*y. Expanding
P@x - alpha*y == 0 implies that alpha**2 == (x.T@x)/(y.T@y)
"""
alpha = np.linalg.norm(x) / np.linalg.norm(y)
self.v = x - alpha*y
denom = self.v.T @ self.v
beta = (0 if denom == 0 else 2.0 / denom)
self.beta_times_v = beta * self.v
def get_matrix(self):
"""
Return the Householder matrix P == I - beta*(v @ v.T)
"""
return np.identity(self.v.size) - np.outer(self.beta_times_v, self.v)
def pre_mult(self, A):
"""
Compute the matrix product P @ A:
P @ A == (I - beta * v @ v.T) @ A == A - (beta * v) @ (v.T @ A)
"""
return A - np.outer(self.beta_times_v, self.v.T @ A).reshape(A.shape)
def __matmul__(self, A): return self.pre_mult(A)
def qr(A, verbose=False):
"""
Compute the QR decomposition, A = Q @ R, where Q @ Q.T = I, and R is upper
triangular, using Householder transformations.
"""
m, n = A.shape
y = np.zeros(m)
y[0] = 1
R = A.copy()
k = min(m, n)
h_list = [None] * k
if verbose: f = open("qr.txt", "w")
# Iterate through columns of A
for i in range(k):
# Apply Householder transform to sub-matrix to remove sub-diagonal
if verbose: print("i={}\n".format(i), R, file=f, end="\n\n")
h = Householder(R[i:, i], y[:m-i])
h_list[i] = h
R[i:, i:] = h @ R[i:, i:]
# Create Q matrix from list of Householder transformations
Q = np.identity(A.shape[0])
for i in reversed(range(k)):
Q[i:, i:] = h_list[i] @ Q[i:, i:]
if verbose:
print("A=", A, "R=", R, "Q=", Q, "Errors=", qr_errors(A, Q, R),
file=f, sep="\n\n")
return Q, R
def qr_errors(A, Q, R):
recon_error = np.max(np.abs(A - Q @ R))
ortho_error = np.max(np.abs(np.identity(A.shape[0]) - Q @ Q.T))
triang_error = np.max(np.abs(
[e for i in range(R.shape[1]) for e in R[i+1:, i]]))
return recon_error, ortho_error, triang_error""" Module to test the Householder class for Householder transformations """
from householder import Householder
import numpy as np
from numpy.linalg import norm
np.random.seed(0)
TOL = 1e-15
N = 10
M = 12
def gaussian_vector(m):
""" Generate an m-dimensional Gaussian vector """
return np.random.normal(size=m)
def gaussian_matrix(m, n, r=None):
"""
Generate an m*n matrix using Gaussian random numbers. If r is None, the
matrix will almost surely be full rank; otherwise, if r is a non-negative
integer, the matrix will have rank r.
"""
if r is None:
return np.random.normal(size=[m, n])
else:
A = np.zeros([m, n])
for _ in range(r):
A += np.outer(gaussian_vector(m), gaussian_vector(n))
return A
def test_symmetric():
""" Test that the Householder matrix is symmetric """
x, y = gaussian_vector(N), gaussian_vector(N)
h = Householder(x, y)
P = h.get_matrix()
assert np.allclose(P, P.T, 0, TOL)
def test_orthogonal():
""" Test that the Householder matrix is orthogonal """
x, y = gaussian_vector(N), gaussian_vector(N)
h = Householder(x, y)
P = h.get_matrix()
assert np.allclose(P @ P.T, np.identity(N), 0, TOL)
def test_self_inverse():
""" Test that the Householder transformation is self-inverse """
x, y = gaussian_vector(N), gaussian_vector(N)
h = Householder(x, y)
assert np.allclose(x, h @ (h @ x), 0, TOL)
P = h.get_matrix()
assert np.allclose(P @ P, np.identity(N), 0, TOL)
def test_householder_multiplication():
"""
Test that the efficient Householder transform implementation is equivalent
to matrix multiplication
"""
x, y = gaussian_vector(N), gaussian_vector(N)
h = Householder(x, y)
z = gaussian_matrix(N, M)
assert np.allclose(h @ z, h.get_matrix() @ z, 0, TOL)
def test_householder_direction():
"""
Test that the Householder transform points vectors in the correct direction
"""
x, y = gaussian_vector(N), gaussian_vector(N)
h = Householder(x, y)
z = h @ x
assert np.allclose(z/norm(z), y/norm(y), 0, TOL)
def test_zero_input():
"""
Verify that if x is a zero-vector, then the Householder transform for any
given y is the identity transformation
"""
x, y = np.zeros(N), gaussian_vector(N)
h = Householder(x, y)
assert np.allclose(h.get_matrix(), np.identity(N), 0, TOL)
def test_zero_out_column():
"""
Test using a Householder transform to zero out the sub-diagonal of the first
column of a matrix
"""
A = gaussian_matrix(M, N)
y = np.zeros(M)
y[0] = 1
h = Householder(A[:, 0], y)
B = h @ A
assert np.allclose(B[1:, 0], 0, 0, TOL)""" Module to test the QR decomposition using Householder transformations """
from householder import qr, qr_errors
from test_householder import gaussian_matrix
import numpy as np
np.random.seed(0)
np.set_printoptions(precision=10, linewidth=1000, suppress=True,
threshold=1000)
TOL = 1e-13
N = 13
M = 17
def _test_qr(m, n):
"""
Test that the QR decomposition works for m*n matrices, including both
full-rank and low-rank matrics
"""
A = gaussian_matrix(m, n)
Q, R = qr(A)
for e in qr_errors(A, Q, R): assert e < TOL
for rank in range(4):
A = gaussian_matrix(m, n, rank)
Q, R = qr(A)
for e in qr_errors(A, Q, R): assert e < TOL
def test_qr_square():
"""
Test that the QR decomposition works for square matrices, including both
full-rank and low-rank matrics
"""
_test_qr(M, M)
def test_qr_skinny():
"""
Test that the QR decomposition works for skinny matrices (more rows than
columns), including both full-rank and low-rank matrics
"""
_test_qr(M, N)
def test_qr_fat():
"""
Test that the QR decomposition works for fat matrices (more columns than
rows), including both full-rank and low-rank matrics
"""
_test_qr(N, M)
def test_small_example():
""" Test small 3x3 example from the Wikipedia page on QR decomposition """
A = np.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]])
Q, R = qr(A)
print("A=", A, "Q=", Q, "R=", R, sep="\n\n",
file=open("Small QR example.txt", "w"))
for e in qr_errors(A, Q, R): assert e < TOL
def test_integer_matrix():
""" Test randomly generated integer matrix, and print output to file """
A = np.random.choice(4, [10, 14]).astype(np.float)
Q, R = qr(A, verbose=True)
for e in qr_errors(A, Q, R): assert e < TOL============================= test session starts =============================
platform win32 -- Python 3.7.6, pytest-5.4.1, py-1.8.1, pluggy-0.13.1
rootdir: c:\Users\Jake\Documents\Programming\Linear algebra
collected 12 items
test_householder.py ....... [ 58%]
test_qr.py ..... [100%]
- generated xml file: C:\Users\Jake\AppData\Local\Temp\tmp-16828fDgc7KQLfgsp.xml -
============================= 12 passed in 0.18s ==============================
After running qr_test.py with pytest, the content of Small QR example.txt is:
A=
[[ 12 -51 4]
[ 6 167 -68]
[ -4 24 -41]]
Q=
[[ 0.8571428571 -0.3942857143 -0.3314285714]
[ 0.4285714286 0.9028571429 0.0342857143]
[-0.2857142857 0.1714285714 -0.9428571429]]
R=
[[ 14 21 -14]
[ 0 175 -70]
[ 0 0 35]]
After running qr_test.py with pytest, the content of qr.txt is:
i=0
[[1. 1. 2. 1. 2. 1. 1. 2. 0. 1. 1. 3. 2. 2.]
[0. 3. 0. 0. 1. 1. 2. 3. 2. 3. 2. 0. 0. 0.]
[1. 3. 3. 1. 3. 2. 0. 1. 0. 3. 1. 3. 0. 3.]
[1. 3. 1. 1. 1. 2. 3. 3. 0. 1. 1. 0. 3. 3.]
[3. 1. 1. 0. 2. 2. 1. 3. 3. 2. 3. 2. 1. 1.]
[3. 3. 2. 3. 1. 3. 3. 3. 1. 3. 2. 0. 2. 2.]
[2. 1. 2. 3. 1. 0. 2. 0. 3. 3. 2. 3. 3. 3.]
[1. 0. 2. 0. 1. 3. 3. 1. 3. 2. 2. 0. 1. 1.]
[2. 3. 1. 1. 3. 1. 3. 0. 0. 3. 0. 3. 3. 1.]
[3. 3. 3. 3. 3. 2. 2. 2. 0. 1. 0. 0. 2. 2.]]
i=1
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 3. 0. 0. 1. 1. 2. 3. 2. 3. 2. 0. 0. 0. ]
[ 0. 2.0915890648 2.4043656599 0.3052968826 2.3738359716 1.2442375061 -0.877881247 0.4348953481 -0.6411234534 2.0610593765 0.4579453239 2.839261008 -0.6261640284 2.4043656599]
[ 0. 2.0915890648 0.4043656599 0.3052968826 0.3738359716 1.2442375061 2.122118753 2.4348953481 -0.6411234534 0.0610593765 0.4579453239 -0.160738992 2.3738359716 2.4043656599]
[ 0. -1.7252328057 -0.7869030204 -2.0841093523 0.1215079148 -0.2672874818 -1.6336437409 1.3046860444 1.0766296397 -0.8168218705 1.3738359716 1.517783024 -0.8784920852 -0.7869030204]
[ 0. 0.2747671943 0.2130969796 0.9158906477 -0.8784920852 0.7327125182 0.3663562591 1.3046860444 -0.9233703603 0.1831781295 0.3738359716 -0.482216976 0.1215079148 0.2130969796]
[ 0. -0.8168218705 0.8087313197 1.6105937652 -0.2523280568 -1.5115249879 0.2442375061 -1.1302093037 1.7177530932 1.122118753 0.9158906477 2.678522016 1.7476719432 1.8087313197]
[ 0. -0.9084109352 1.4043656599 -0.6947031174 0.3738359716 2.2442375061 2.122118753 0.4348953481 2.3588765466 1.0610593765 1.4579453239 -0.160738992 0.3738359716 0.4043656599]
[ 0. 1.1831781295 -0.1912686803 -0.3894062348 1.7476719432 -0.5115249879 1.2442375061 -1.1302093037 -1.2822469068 1.122118753 -1.0841093523 2.678522016 1.7476719432 -0.1912686803]
[ 0. 0.2747671943 1.2130969796 0.9158906477 1.1215079148 -0.2672874818 -0.6336437409 0.3046860444 -1.9233703603 -1.8168218705 -1.6261640284 -0.482216976 0.1215079148 0.2130969796]]
i=2
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 1.1549789687 -0.6625378782 1.0783417646 0.6702201207 -1.0932752482 0.9476616396 2.5924936929 2.2394827507 2.2662058224 1.9594112255 -1.8931477152 0.2399351948]
[ 0. 0. -0.8450210313 -0.6625378782 -0.9216582354 0.6702201207 1.9067247518 2.9476616396 2.5924936929 0.2394827507 2.2662058224 -1.0405887745 1.1068522848 0.2399351948]
[ 0. -0. 0.2436449696 -1.285797529 1.1900873711 0.2061867971 -1.4559774854 0.8817343075 -1.5905971476 -0.9639931593 -0.1176952215 2.2435210451 0.1665706653 0.9984125116]
[ 0. 0. 0.0489679537 0.7887484446 -1.0486781483 0.6573051884 0.3380604498 1.3720469604 -0.4985775911 0.2066171935 0.6113829507 -0.5978007971 -0.0449327875 -0.0712392202]
[ 0. -0. 1.2966504771 1.9885593096 0.2535973275 -1.2873557367 0.3283546656 -1.3304583529 0.4549384943 1.0524396193 0.2097162143 3.0221270564 2.2424632014 2.6539997191]
[ 0. -0. 1.9469944925 -0.2743568386 0.9364900436 2.4935425338 2.215667849 0.2121926604 0.954464358 0.9835672214 0.6725885641 0.2213939887 0.9241074639 1.3444127925]
[ 0. 0. -0.8980265388 -0.9368947168 1.0148318082 -0.8362373455 1.1223926008 -0.8401457 0.5469580509 1.2230499721 -0.0612056134 2.1808052142 1.0309597486 -1.4156520127]
[ 0. 0. 1.0489679537 0.7887484446 0.9513218517 -0.3426948116 -0.6619395502 0.3720469604 -1.4985775911 -1.7933828065 -1.3886170493 -0.5978007971 -0.0449327875 -0.0712392202]]
i=3
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 0.0888696592 -0.7535451914 0.9270768033 2.361645751 2.3491565884 1.5653542522 -0.3929583391 1.3974846306 -1.1277509585 2.0304022029 1.2131258711]
[ 0. -0. 0. -1.5024509138 1.141615326 0.1321272929 -1.5871448871 1.0543013083 -1.2944419651 -0.7816413782 0.1327832169 2.2686525234 -0.0997165435 0.7178123798]
[ 0. 0. 0. 0.7452052804 -1.0584200972 0.6424206521 0.3116983243 1.4067296092 -0.4390560943 0.243266395 0.6617243002 -0.5927498532 -0.0984513959 -0.1276344509]
[ 0. -0. 0. 0.8355549396 -0.0043653157 -1.6814918916 -0.3697031446 -0.4120766423 2.0310422792 2.0228948036 1.5427336259 3.1558738902 0.8253133345 1.1606781236]
[ 0. -0. 0. -2.0056585254 0.5491444524 1.9017246644 1.1674943158 1.5911950701 3.3210738978 2.4407609159 2.6741900786 0.4222224721 -1.2038237208 -0.8978945782]
[ 0. 0. -0. -0.1383537515 1.1934900539 -0.5632688523 1.6058493613 -1.4761930878 -0.5446106114 0.550937845 -0.9844189492 2.088175621 2.012441044 -0.3814161508]
[ 0. 0. 0. -0.1440122769 0.7426345063 -0.6615441879 -1.2266563229 1.1150019823 -0.2235367136 -1.0083012391 -0.3102288706 -0.4896019069 -1.1913826768 -1.279310729 ]]
i=4
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 2.7534136901 -1.4338726432 -1.7009353972 0.0474833282 -1.3870086628 -1.1257264879 -0.6593163351 -1.2623546037 -0.8639488348 1.1818393949 0.7052719388]
[ 0. -0. 0. 0. 0.7580004974 -1.3497244935 -2.8920268505 -1.0524024393 -2.811856078 -0.9318321087 -1.3670148129 2.4174020795 -0.5781938933 0.4314498368]
[ 0. 0. 0. 0. -0.8681497907 1.3774089076 0.9589107693 2.4516401234 0.3135708296 0.3177599604 1.4056137685 -0.6665286061 0.1388700663 0.0143993934]
[ 0. -0. 0. -0. 0.2089736107 -0.8573927023 0.3559781806 0.7595201843 2.874918672 2.106420066 2.3768132227 3.0731501052 1.09140796 1.3199323363]
[ 0. -0. 0. 0. 0.0370476215 -0.0764355813 -0.5744245786 -1.2210953637 1.2954406313 2.240267631 0.672072953 0.620791898 -1.8425548547 -1.280166951 ]
[ 0. 0. -0. 0. 1.1581647397 -0.6997257255 1.48568882 -1.6701896864 -0.6843422547 0.5371074757 -1.1225284089 2.1018732792 1.9683802792 -0.4077859521]
[ 0. 0. 0. 0. 0.7058644232 -0.8035820062 -1.3517313061 0.9130711233 -0.3689832367 -1.0226972562 -0.4539868644 -0.4753440286 -1.2372454813 -1.3067590279]]
i=5
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 2.7534136901 -1.4338726432 -1.7009353972 0.0474833282 -1.3870086628 -1.1257264879 -0.6593163351 -1.2623546037 -0.8639488348 1.1818393949 0.7052719388]
[ 0. -0. 0. 0. 1.7924512427 -2.1080132164 -1.2301554286 -2.288755041 -1.5664988637 -0.3117721178 -1.8719712873 2.8871255981 0.5620029997 -0.4751793256]
[ 0. 0. 0. 0. 0. 0.7410246243 2.3536154796 1.4140467366 1.3587212244 0.8381375171 0.9818353919 -0.2723190449 1.0957659827 -0.7464777397]
[ 0. -0. 0. -0. 0. -0.7042077072 0.0202567621 1.0092808118 2.623338986 1.9811592203 2.4788215379 2.978259343 0.861072139 1.5030841881]
[ 0. -0. 0. 0. 0. -0.0492783738 -0.6339425244 -1.1768168677 1.2508396458 2.2180609209 0.6901573677 0.6039693096 -1.8833896486 -1.2476971072]
[ 0. 0. -0. 0. 0. 0.1492496998 -0.3749322791 -0.2859768609 -2.0786365895 -0.157107905 -0.5571822227 1.5759736952 0.6918227982 0.6072705235]
[ 0. 0. 0. 0. 0. -0.2861585962 -2.4857204477 1.7567046361 -1.2187611074 -1.4457993612 -0.1094264181 -0.7958630306 -2.0152664996 -0.6881145215]]
i=6
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 2.7534136901 -1.4338726432 -1.7009353972 0.0474833282 -1.3870086628 -1.1257264879 -0.6593163351 -1.2623546037 -0.8639488348 1.1818393949 0.7052719388]
[ 0. -0. 0. 0. 1.7924512427 -2.1080132164 -1.2301554286 -2.288755041 -1.5664988637 -0.3117721178 -1.8719712873 2.8871255981 0.5620029997 -0.4751793256]
[ 0. 0. 0. 0. 0. 1.0731339906 2.2517350862 -0.1400443446 -0.8047885149 -0.4594858469 -1.0286663557 -1.738752378 0.9116909015 -1.1765656053]
[ 0. -0. 0. -0. 0. 0. -0.1957713464 -2.2860279871 -1.9641866529 -0.770333109 -1.7842645701 -0.1311791731 0.4707576704 0.5911219842]
[ 0. -0. 0. 0. 0. 0. -0.6490595329 -1.4074128367 0.9298181619 2.0255196209 0.3918392032 0.3863799935 -1.9107027013 -1.3115135266]
[ 0. 0. -0. 0. 0. -0. -0.3291473068 0.4124305038 -1.1063569376 0.4260430602 0.3463357526 2.2349877337 0.7745460014 0.8005516873]
[ 0. 0. 0. 0. 0. 0. -2.5735046336 0.4176381501 -3.0829268691 -2.5638830947 -1.8417544282 -2.0594001258 -2.173872887 -1.0586952769]]
i=7
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 2.7534136901 -1.4338726432 -1.7009353972 0.0474833282 -1.3870086628 -1.1257264879 -0.6593163351 -1.2623546037 -0.8639488348 1.1818393949 0.7052719388]
[ 0. -0. 0. 0. 1.7924512427 -2.1080132164 -1.2301554286 -2.288755041 -1.5664988637 -0.3117721178 -1.8719712873 2.8871255981 0.5620029997 -0.4751793256]
[ 0. 0. 0. 0. 0. 1.0731339906 2.2517350862 -0.1400443446 -0.8047885149 -0.4594858469 -1.0286663557 -1.738752378 0.9116909015 -1.1765656053]
[ 0. -0. 0. -0. 0. 0. 2.6815795244 0.0561196825 3.0128157823 1.9742321199 1.7604359441 1.6181258242 2.4192947096 1.1920524265]
[ 0. -0. 0. 0. 0. 0. 0. -0.8790819725 2.0525074066 2.6446259813 1.1914363894 0.7809801432 -1.4711607069 -1.1759584102]
[ 0. 0. -0. 0. 0. -0. -0. 0.6803545882 -0.5370252777 0.7400006879 0.7518228413 2.4350950736 0.997443989 0.8692936119]
[ 0. 0. 0. 0. 0. 0. 0. 2.5124566656 1.3685075937 -0.1091420323 1.3286281936 -0.4948205114 -0.4310999948 -0.5212226685]]
i=8
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 2.7534136901 -1.4338726432 -1.7009353972 0.0474833282 -1.3870086628 -1.1257264879 -0.6593163351 -1.2623546037 -0.8639488348 1.1818393949 0.7052719388]
[ 0. -0. 0. 0. 1.7924512427 -2.1080132164 -1.2301554286 -2.288755041 -1.5664988637 -0.3117721178 -1.8719712873 2.8871255981 0.5620029997 -0.4751793256]
[ 0. 0. 0. 0. 0. 1.0731339906 2.2517350862 -0.1400443446 -0.8047885149 -0.4594858469 -1.0286663557 -1.738752378 0.9116909015 -1.1765656053]
[ 0. -0. 0. -0. 0. 0. 2.6815795244 0.0561196825 3.0128157823 1.9742321199 1.7604359441 1.6181258242 2.4192947096 1.1920524265]
[ 0. -0. 0. 0. 0. 0. 0. 2.7473816584 0.4617582529 -0.7627607111 1.0199735522 -0.0993791853 0.3234957882 0.1148884311]
[ 0. 0. -0. 0. 0. -0. -0. 0. -0.2385875721 1.3792546294 0.7839906885 2.6002577691 0.6607515787 0.627120062 ]
[ 0. 0. 0. 0. 0. 0. 0. 0. 2.4705973212 2.2515355532 1.4474196643 0.115102788 -1.6744591271 -1.4155365458]]
i=9
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 2.7534136901 -1.4338726432 -1.7009353972 0.0474833282 -1.3870086628 -1.1257264879 -0.6593163351 -1.2623546037 -0.8639488348 1.1818393949 0.7052719388]
[ 0. -0. 0. 0. 1.7924512427 -2.1080132164 -1.2301554286 -2.288755041 -1.5664988637 -0.3117721178 -1.8719712873 2.8871255981 0.5620029997 -0.4751793256]
[ 0. 0. 0. 0. 0. 1.0731339906 2.2517350862 -0.1400443446 -0.8047885149 -0.4594858469 -1.0286663557 -1.738752378 0.9116909015 -1.1765656053]
[ 0. -0. 0. -0. 0. 0. 2.6815795244 0.0561196825 3.0128157823 1.9742321199 1.7604359441 1.6181258242 2.4192947096 1.1920524265]
[ 0. -0. 0. 0. 0. 0. 0. 2.7473816584 0.4617582529 -0.7627607111 1.0199735522 -0.0993791853 0.3234957882 0.1148884311]
[ 0. 0. -0. 0. 0. -0. -0. 0. 2.4820908833 2.1085306457 1.3653572208 -0.135376408 -1.7302192187 -1.4692628202]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 1.5892936155 0.9194915681 2.5992811209 0.4967368241 0.4881496184]]
A=
[[1. 1. 2. 1. 2. 1. 1. 2. 0. 1. 1. 3. 2. 2.]
[0. 3. 0. 0. 1. 1. 2. 3. 2. 3. 2. 0. 0. 0.]
[1. 3. 3. 1. 3. 2. 0. 1. 0. 3. 1. 3. 0. 3.]
[1. 3. 1. 1. 1. 2. 3. 3. 0. 1. 1. 0. 3. 3.]
[3. 1. 1. 0. 2. 2. 1. 3. 3. 2. 3. 2. 1. 1.]
[3. 3. 2. 3. 1. 3. 3. 3. 1. 3. 2. 0. 2. 2.]
[2. 1. 2. 3. 1. 0. 2. 0. 3. 3. 2. 3. 3. 3.]
[1. 0. 2. 0. 1. 3. 3. 1. 3. 2. 2. 0. 1. 1.]
[2. 3. 1. 1. 3. 1. 3. 0. 0. 3. 0. 3. 3. 1.]
[3. 3. 3. 3. 3. 2. 2. 2. 0. 1. 0. 0. 2. 2.]]
R=
[[ 6.2449979984 5.764613537 5.1241009218 4.6437164603 5.2842290756 4.963972768 5.6044853832 4.963972768 3.3626912299 5.9247416908 3.8430756913 3.8430756913 5.2842290756 5.1241009218]
[ 0. 4.8753698085 1.1202306032 0.8677842701 2.1615717274 1.5146779987 2.193127519 2.5402412271 -0.8993400618 2.8400212476 0.3786694997 0.788894791 1.1360084991 1.9406811859]
[ 0. 0. 3.080368994 1.049552737 1.4613892465 1.2554709918 -0.0567325804 -0.4160389228 0.2521448017 0.7984585387 0.2868147119 1.7608111985 0.2111712714 2.4573612131]
[ 0. 0. -0. 2.7534136901 -1.4338726432 -1.7009353972 0.0474833282 -1.3870086628 -1.1257264879 -0.6593163351 -1.2623546037 -0.8639488348 1.1818393949 0.7052719388]
[ 0. -0. 0. 0. 1.7924512427 -2.1080132164 -1.2301554286 -2.288755041 -1.5664988637 -0.3117721178 -1.8719712873 2.8871255981 0.5620029997 -0.4751793256]
[ 0. 0. 0. 0. 0. 1.0731339906 2.2517350862 -0.1400443446 -0.8047885149 -0.4594858469 -1.0286663557 -1.738752378 0.9116909015 -1.1765656053]
[ 0. -0. 0. -0. 0. 0. 2.6815795244 0.0561196825 3.0128157823 1.9742321199 1.7604359441 1.6181258242 2.4192947096 1.1920524265]
[ 0. -0. 0. 0. 0. 0. 0. 2.7473816584 0.4617582529 -0.7627607111 1.0199735522 -0.0993791853 0.3234957882 0.1148884311]
[ 0. 0. -0. 0. 0. -0. -0. 0. 2.4820908833 2.1085306457 1.3653572208 -0.135376408 -1.7302192187 -1.4692628202]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 1.5892936155 0.9194915681 2.5992811209 0.4967368241 0.4881496184]]
Q=
[[ 0.1601281538 0.0157778958 0.3771665992 -0.0556174288 0.2727023011 0.1751565529 0.0123284654 0.6889502408 -0.1889948857 0.4617274469]
[ 0. 0.615337937 -0.2237785115 -0.108633597 -0.0886140254 -0.0211247362 0.2168522206 0.3549463231 0.6101634059 -0.0766242973]
[ 0.1601281538 0.4260031871 0.552617357 -0.2517855092 0.0359241527 -0.4533141359 -0.2697934352 -0.3503200136 0.0354270647 0.1433402804]
[ 0.1601281538 0.4260031871 -0.0966555073 -0.0042941445 -0.3525318923 -0.064509386 0.3261488956 0.0883066916 -0.7104344049 -0.1902396348]
[ 0.4803844614 -0.3628916039 -0.342496689 -0.5652580553 -0.0357233656 -0.411625728 -0.0022761639 0.1715913416 0.0196112464 -0.0132110857]
[ 0.4803844614 0.0473336875 -0.1670459311 0.3281304431 -0.5167004177 0.2071763816 -0.3443128649 -0.0922683192 0.1111840056 0.4267180697]
[ 0.3202563076 -0.173556854 0.1796531712 0.5356532324 0.105086931 -0.3911664406 0.5894266982 -0.0649745083 0.1727057896 0.0508626801]
[ 0.1601281538 -0.1893347498 0.4517594364 -0.382591756 -0.3602203638 0.4795536939 0.387324201 -0.1585721914 0.1912096526 -0.128808086 ]
[ 0.3202563076 0.2366684373 -0.2941689353 -0.1393945358 0.5724773916 0.3641570496 0.2089309888 -0.4211763418 -0.0773376022 0.2193040234]
[ 0.4803844614 0.0473336875 0.1575905011 0.2043847607 0.2354228357 0.1768266988 -0.3376521382 0.1553226351 0.0298901474 -0.6889581218]]
Errors=
(1.3322676295501878e-15, 2.748003089305132e-16, 4.440892098500626e-16)