import numpy as np
import scipy
from scipy import linalg
#https://gist.github.com/samubernard/746c684771bc74d446ec
There are a several different algorithms for calculating the matrices $Q$ and $R$.
We will outline the method of Householder Reflections, which is known to be more numerically stable than the alternative Gramm-Schmidt method.
I'm sure I got some of the code from this repo
def QR(A):
N = A.shape[0]
Q = A.astype(float).T
for i in range(N):
Q[i] /= np.linalg.norm(Q[i])
for j in range(i+1,N):
Q[j] = Q[j] - (Q[i] @ Q[j]) * Q[i]
return Q.T, (Q@A) * np.tri(N).T
QR(A)
(array([[ 0.85714286, -0.39428571, -0.33142857],
[ 0.42857143, 0.90285714, 0.03428571],
[-0.28571429, 0.17142857, -0.94285714]]),
array([[ 14., 21., -14.],
[ -0., 175., -70.],
[ -0., -0., 35.]]))
def MGS(A):
"""factorisation QR Gram-Schmidt modifiée - stable
Note:
Cette factorisation est stable numériquement
Tire de: LN Trefethen & D Bau III, Numerical Linear Algebra, 1997 SIAM Philadelphia
"""
n = A.shape[0]
R = np.zeros_like(A)
Q = np.zeros_like(A)
V = A.copy()
for i in range(n):
# diagonal elements are norms
R[i,i] = np.linalg.norm(V[:,i])
Q[:,i] = V[:,i]/R[i,i]
for j in range(i+1, n):
R[i,j] = Q[:,i] @ V[:,j]
#make ortho to the other columns
V[:,j] = V[:,j] - R[i,j]*Q[:,i]
return Q,R
A = np.array([[2,1,1],[1,1,1],[0,1,1]]).astype(float)
B = np.array([[12, -51, 4],
[6, 167, -68],
[-4, 24, -41]]).T
Q,R = MGS(A)
Q @ R
array([[2., 1., 1.],
[1., 1., 1.],
[0., 1., 1.]])
You may have to install tabulate in the terminal or with conda
! pip install tabulate
# A is a square random matrix of size n
n = 5
A = np.random.rand(n, n)
#make symmetric so eigenvalues real and in fact >= 0
A = A.T@A
np.max(np.abs(np.tril(A, k = -1)))
1.8319623963241725
np.linalg.norm(np.tril(A, k = -1).ravel(), ord=np.inf)
1.8319623963241725
? np.linalg.norm
Signature: np.linalg.norm(x, ord=None, axis=None, keepdims=False) Docstring: Matrix or vector norm. This function is able to return one of eight different matrix norms, or one of an infinite number of vector norms (described below), depending on the value of the ``ord`` parameter. Parameters ---------- x : array_like Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord` is None. If both `axis` and `ord` are None, the 2-norm of ``x.ravel`` will be returned. ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional Order of the norm (see table under ``Notes``). inf means numpy's `inf` object. The default is None. axis : {None, int, 2-tuple of ints}, optional. If `axis` is an integer, it specifies the axis of `x` along which to compute the vector norms. If `axis` is a 2-tuple, it specifies the axes that hold 2-D matrices, and the matrix norms of these matrices are computed. If `axis` is None then either a vector norm (when `x` is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default is None. .. versionadded:: 1.8.0 keepdims : bool, optional If this is set to True, the axes which are normed over are left in the result as dimensions with size one. With this option the result will broadcast correctly against the original `x`. .. versionadded:: 1.10.0 Returns ------- n : float or ndarray Norm of the matrix or vector(s). See Also -------- scipy.linalg.norm : Similar function in SciPy. Notes ----- For values of ``ord < 1``, the result is, strictly speaking, not a mathematical 'norm', but it may still be useful for various numerical purposes. The following norms can be calculated: ===== ============================ ========================== ord norm for matrices norm for vectors ===== ============================ ========================== None Frobenius norm 2-norm 'fro' Frobenius norm -- 'nuc' nuclear norm -- inf max(sum(abs(x), axis=1)) max(abs(x)) -inf min(sum(abs(x), axis=1)) min(abs(x)) 0 -- sum(x != 0) 1 max(sum(abs(x), axis=0)) as below -1 min(sum(abs(x), axis=0)) as below 2 2-norm (largest sing. value) as below -2 smallest singular value as below other -- sum(abs(x)**ord)**(1./ord) ===== ============================ ========================== The Frobenius norm is given by [1]_: :math:`||A||_F = [\sum_{i,j} abs(a_{i,j})^2]^{1/2}` The nuclear norm is the sum of the singular values. Both the Frobenius and nuclear norm orders are only defined for matrices and raise a ValueError when ``x.ndim != 2``. References ---------- .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*, Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15 Examples -------- >>> from numpy import linalg as LA >>> a = np.arange(9) - 4 >>> a array([-4, -3, -2, ..., 2, 3, 4]) >>> b = a.reshape((3, 3)) >>> b array([[-4, -3, -2], [-1, 0, 1], [ 2, 3, 4]]) >>> LA.norm(a) 7.745966692414834 >>> LA.norm(b) 7.745966692414834 >>> LA.norm(b, 'fro') 7.745966692414834 >>> LA.norm(a, np.inf) 4.0 >>> LA.norm(b, np.inf) 9.0 >>> LA.norm(a, -np.inf) 0.0 >>> LA.norm(b, -np.inf) 2.0 >>> LA.norm(a, 1) 20.0 >>> LA.norm(b, 1) 7.0 >>> LA.norm(a, -1) -4.6566128774142013e-010 >>> LA.norm(b, -1) 6.0 >>> LA.norm(a, 2) 7.745966692414834 >>> LA.norm(b, 2) 7.3484692283495345 >>> LA.norm(a, -2) 0.0 >>> LA.norm(b, -2) 1.8570331885190563e-016 # may vary >>> LA.norm(a, 3) 5.8480354764257312 # may vary >>> LA.norm(a, -3) 0.0 Using the `axis` argument to compute vector norms: >>> c = np.array([[ 1, 2, 3], ... [-1, 1, 4]]) >>> LA.norm(c, axis=0) array([ 1.41421356, 2.23606798, 5. ]) >>> LA.norm(c, axis=1) array([ 3.74165739, 4.24264069]) >>> LA.norm(c, ord=1, axis=1) array([ 6., 6.]) Using the `axis` argument to compute matrix norms: >>> m = np.arange(8).reshape(2,2,2) >>> LA.norm(m, axis=(1,2)) array([ 3.74165739, 11.22497216]) >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :]) (3.7416573867739413, 11.224972160321824) File: ~/anaconda3/lib/python3.9/site-packages/numpy/linalg/linalg.py Type: function
import numpy as np
#! pip install tabulate
from tabulate import tabulate
# A is a square random matrix of size n
n = 5
A = np.random.rand(n, n)
#make symmetric so eigenvalues real and in fact >= 0
A = A.T@A
def eigen_qr_simple(A, max_iter=50000):
A_k = np.copy(A)
n = A.shape[0]
QQ = np.identity(n)
peek_time = max_iter // 5
TT = []
for k in range(max_iter):
Q, R = np.linalg.qr(A_k)
A_k = R @ Q
QQ = QQ @ Q
# we "peek" into the structure of matrix A_k from time to time
# to see how it looks
if k % peek_time == 0:
print("A",k,"=")
print(tabulate(A_k))
print("\n")
# break if the lower triangular entries all close to 0
# use the sup norm
TT.append(np.linalg.norm(np.tril(A, k = -1).ravel(), ord=np.inf))
if TT[-1] < .001: break
# return diagonal elts and the errors
return [A_k[k,k] for k in range(A_k.shape[0] )]
print(eigen_qr_simple(A))
# We compare our results with the official numpy algorithm
print(np.linalg.eigvals(A))
A 0 = ------------ ----------- ------------ ------------ ------------ 6.56293 -0.0358214 0.0557912 0.0249674 -9.08991e-05 -0.0358214 0.290227 -0.0630039 -0.0575942 2.4846e-05 0.0557912 -0.0630039 0.212473 0.242525 0.000157138 0.0249674 -0.0575942 0.242525 0.382462 0.000159074 -9.08991e-05 2.4846e-05 0.000157138 0.000159074 0.0001425 ------------ ----------- ------------ ------------ ------------ A 10000 = ------- ------------- ------------ ----------- ------------ 6.56375 3.63751e-16 7.28712e-16 6.53224e-16 7.43012e-16 0 0.57802 -8.53102e-17 1.36493e-16 -2.33506e-16 0 -4.94066e-324 0.267342 1.63719e-16 -1.09348e-16 0 0 0 0.0389769 -2.79282e-16 0 0 0 0 0.000142356 ------- ------------- ------------ ----------- ------------ A 20000 = ------- ------------- ------------ ----------- ------------ 6.56375 3.63751e-16 7.28712e-16 6.53224e-16 7.43012e-16 0 0.57802 -8.53102e-17 1.36493e-16 -2.33506e-16 0 -4.94066e-324 0.267342 1.63719e-16 -1.09348e-16 0 0 0 0.0389769 -2.79282e-16 0 0 0 0 0.000142356 ------- ------------- ------------ ----------- ------------ A 30000 = ------- ------------- ------------ ----------- ------------ 6.56375 3.63751e-16 7.28712e-16 6.53224e-16 7.43012e-16 0 0.57802 -8.53102e-17 1.36493e-16 -2.33506e-16 0 -4.94066e-324 0.267342 1.63719e-16 -1.09348e-16 0 0 0 0.0389769 -2.79282e-16 0 0 0 0 0.000142356 ------- ------------- ------------ ----------- ------------ A 40000 = ------- ------------- ------------ ----------- ------------ 6.56375 3.63751e-16 7.28712e-16 6.53224e-16 7.43012e-16 0 0.57802 -8.53102e-17 1.36493e-16 -2.33506e-16 0 -4.94066e-324 0.267342 1.63719e-16 -1.09348e-16 0 0 0 0.0389769 -2.79282e-16 0 0 0 0 0.000142356 ------- ------------- ------------ ----------- ------------ [6.563748950024317, 0.5780196926901723, 0.2673422259774415, 0.0389769426396412, 0.00014235649279067998] [6.56374895e+00 5.78019693e-01 1.42356493e-04 3.89769426e-02 2.67342226e-01]
def horner(x, P):
val = 0
for coeff in reversed(P):
val *= x
val += coeff
return val
def companion(P):
n = len(P) - 1
C = np.zeros((n,n))
C[1:,:-1] = np.identity(n-1)
C[:,-1] = [-x/P[-1] for x in P[:-1]]
return C
P = [6,-5,-2,1]
[ (x, horner( np.round(x) ,P)) for x in np.linalg.eigvals( companion(P)) ]
[(-2.000000000000001, 0.0), (1.0000000000000016, 0.0), (2.999999999999999, 0.0)]