Demo of Matrix Norms

In [1]:
import numpy as np
from numpy import cos, sin
import matplotlib.pyplot as plt
%matplotlib inline

Generate a random matrix to play with

In [2]:
# Generate a random matrix to play with...
A = np.random.rand(2,2)-0.5
A *= 2.

# ... Or use the one from the lectures.
A = np.array([[1, 2],[0.25, -1]], dtype=float)
print(A)

[[ 1.    2.  ]
 [ 0.25 -1.  ]]

Apply it to a bunch of unit vectors

In [3]:
S = []
X = []
thetas = np.linspace(0., 2.*np.pi, 360)
for theta in thetas:
    s = np.array([np.cos(theta), np.sin(theta)])
    S.append( s )
    X.append( A @ s )
S = np.array(S)
X = np.array(X)
In [4]:
plt.plot(S[:,0], S[:,1], 'b')
plt.plot(X[:,0], X[:,1], 'r')
plt.axis('equal');

What is the longest mapped vector?

In [5]:
normX = [np.linalg.norm(x) for x in X]
plt.plot(thetas, normX)
normA = np.max(normX)
plt.xlabel(r'$\theta$'); plt.ylabel(r'$\||Ax\||_2$')
plt.plot([0,thetas[-1]], [normA, normA], 'r--')
plt.ylim([0, normA*1.1])
plt.title(r'$\||A\||_2 = $'+str(normA));

Condition Number

In [9]:
A1 = np.array([[1,1],[1,-1]])
A2 = np.array([[1,1],[1,0.9]])
b = np.array([0,0])
x1_true = np.linalg.solve(A1, b)
x2_true = np.linalg.solve(A2, b)
#A1[0,0]*x[0] + A1[0,1]*x[1] = b[0]
In [10]:
x = np.array([-10, 10])
y1p = ( b[0] - A1[0,1]*x )/A1[0,0]
y1m = ( b[1] - A1[1,1]*x )/A1[1,0]
y2p = ( b[0] - A2[0,1]*x )/A2[0,0]
y2m = ( b[1] - A2[1,1]*x )/A2[1,0]
In [12]:
plt.figure(figsize=[13,6])
trials = 500
sig = 0.5

maxK2 = 0.
maxK1 = 0.
X1 = []
X2 = []
for idx in range(trials):
    db = np.random.normal(scale=sig, size=(2,))
    x1 = np.linalg.solve(A1,b+db)
    x2 = np.linalg.solve(A2,b+db)
    dx1 = x1 - x1_true
    dx2 = x2 - x2_true
    maxK1 = max(maxK1, np.linalg.norm(dx1) / np.linalg.norm(db))
    maxK2 = max(maxK2, np.linalg.norm(dx2) / np.linalg.norm(db))
    X1.append(x1)
    X2.append(x2)
X1 = np.array(X1); X2 = np.array(X2)
plt.subplot(1,2,1); plt.plot(X1[:,0], X1[:,1], 'b.');
plt.axis(10.*np.array([-1,1,-1,1])); plt.grid(True);
plt.subplot(1,2,2); plt.plot(X2[:,0], X2[:,1], 'r.');
plt.axis(10.*np.array([-1,1,-1,1])); plt.grid(True);
plt.subplot(1,2,1); plt.plot(x, y1p, 'k'); plt.plot(x, y1m, 'k')
plt.subplot(1,2,2); plt.plot(x, y2p, 'k'); plt.plot(x, y2m, 'k');
In [13]:
print(maxK1)
print(np.linalg.cond(A1))

0.7071067811865478
1.0000000000000004
In [14]:
print(maxK2)
print(np.linalg.cond(A2))

19.511965080466783
38.073735174775784
In [ ]: