1import numpy as np
2import Tools
3from numpy.linalg import qr
4import math
5
6def randomIsometry(rows,cols,rank):
7    if rank==1:
8       r=np.random.randn(rows)
9       r = Tools.normalize(r)[np.newaxis]
10       c=np.random.randn(cols)
11       c = Tools.normalize(c)[np.newaxis]
12       result=r.T.dot(c)
13    else:
14        a = np.random.randn(rows,rows)
15        b = np.random.randn(cols,cols)
16
17        diagDim = min(rows,cols)
18        d = np.zeros((rows,cols))
19
20        diag = np.ones(diagDim)
21        diag[rank:] = 0
22        np.fill_diagonal(d,diag)
23
24        qa,_ = qr(a)
25        qb,_ = qr(b)
26
27        result = qa .dot(d.dot(qb))
28    return(result)
29
30def kahan_matrix(rows):
31    cols = rows
32    eps = np.finfo(np.float32).eps
33    s = math.pow(eps,1.0/rows)
34    c = math.sqrt(1-s*s)
35    m = np.zeros((rows,cols))
36    sc = 1.0
37
38    for i in range(rows-1):
39        m[i,i] = sc
40        m[i,i+1:] = - sc * c * np.ones(rows-i-1)
41        sc = sc * s
42
43    m = m + m.T
44
45    return(m)
46
47def householder(x,eps=1e-16):
48    #print(x)
49    v=np.hstack([[1],x[1:]])
50
51    alpha = x[0]
52    xnorm2=x[1:].dot(x[1:])
53    epsilon=eps
54
55    #print(sigma)
56
57    if xnorm2<=epsilon:
58       tau = 0.0
59       v = np.zeros(len(x))
60    else:
61       if np.sign(alpha) <= 0:
62          beta = math.sqrt(alpha*alpha + xnorm2)
63       else:
64          beta = -math.sqrt(alpha*alpha + xnorm2)
65
66       r =  (alpha - beta)
67       v = x / r
68       tau = (beta - alpha) / beta
69       v[0] = 1
70
71    return(v,tau)
72
73def QR(oldm,eps=1e-16):
74    m=np.copy(oldm)
75    (rows,cols) = m.shape
76    hvects=[]
77    tau=[]
78    h=[]
79    for c in range(cols):
80        currentSize = rows - c
81        v,beta=householder(m[c:,c],eps=eps)
82        tau.append(beta)
83        h.append(v)
84        hvects.append((v,beta))
85        t = np.identity(currentSize) - beta * np.outer(v,v)
86        m[c:,c:] = np.dot(t,m[c:,c:])
87
88    hvects.reverse()
89    q=np.identity(rows)
90    c=cols - 1
91    for (v,beta) in hvects:
92       t = np.identity(len(v)) - beta * np.outer(v,v)
93       q[c:,c:] = np.dot(t,q[c:,c:])
94       c = c - 1
95
96    return(q,m,tau,h)