Demo of Matrix Norms¶

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

Generate a random matrix to play with¶

In [13]:
# Generate a random matrix to play with...
A = np.random.rand(2,2)-0.5
A += A.T
#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 [14]:
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 [15]:
k = 25
plt.figure(1)
plt.plot(S[:,0], S[:,1], 'b')
plt.plot(X[:,0], X[:,1], 'r')
plt.plot([0,S[k,0]], [0,S[k,1]], color='lightblue', marker='o')
plt.plot([0,X[k,0]], [0,X[k,1]], color='pink', marker='o')
plt.axis('equal');
In [16]:
from ipywidgets import interact, interactive, fixed, interact_manual
fig = plt.figure(1, figsize=(6,6))
plt.clf()
ax = plt.gca()
plt.plot(S[:,0], S[:,1], 'b')
plt.plot(X[:,0], X[:,1], 'r')
plt.axis('equal');
ax = plt.gca()
unit_line, = ax.plot([0,S[0,0]], [0,S[0,1]], color='lightblue', marker='o')
mapped_line, = ax.plot([0,X[0,0]], [0,X[0,1]], color='pink', marker='o')

def update(theta=0.):
    s = [np.cos(theta), np.sin(theta)]
    unit_line.set_xdata([0, s[0]])
    unit_line.set_ydata([0, s[1]])
    x = A @ s
    mapped_line.set_xdata([0, x[0]])
    mapped_line.set_ydata([0, x[1]])
    plt.title(f'Stretch = {np.linalg.norm(x)}')
    fig.canvas.draw_idle()
    
interact(update, theta=(0, 2.*np.pi, 0.05));

interactive(children=(FloatSlider(value=0.0, description='theta', max=6.283185307179586, step=0.05), Output())…

What is the longest mapped vector?¶

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

Condition Number¶

In [18]:
A1 = np.array([[ 1,1],
               [ 1,-1]])
A2 = np.array([[1,1],
               [1,0.8]])
b = np.array([1,1])
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 [19]:
# Some code to draw the intersecting lines
x = np.array([-10, 10])
y1p = ( b[0] - A1[0,0]*x )/A1[0,1]
y1m = ( b[1] - A1[1,0]*x )/A1[1,1]
y2p = ( b[0] - A2[0,0]*x )/A2[0,1]
y2m = ( b[1] - A2[1,0]*x )/A2[1,1]
In [20]:
plt.figure(figsize=[13,6])
trials = 1000
sig = 0.25

maxK1 = 0.
maxK2 = 0.
X1 = []
X2 = []
for idx in range(trials):
    
    # ADD NOISE TO b, and solve for x =====
    db = np.random.normal(scale=sig, size=(2,))
    relerr_b = np.linalg.norm(db) / np.linalg.norm(b)
    x1 = np.linalg.solve(A1,b+db)
    x2 = np.linalg.solve(A2,b+db)
    # =====
    
    # Relative perturbation in x
    relerr_x1 = np.linalg.norm(x1 - x1_true) / np.linalg.norm(x1_true)
    relerr_x2 = np.linalg.norm(x2 - x2_true) / np.linalg.norm(x2_true)
    
    # Keep track of the maximum displacement
    maxK1 = max(maxK1, relerr_x1 / relerr_b)
    maxK2 = max(maxK2, relerr_x2 / relerr_b)
    X1.append(x1)
    X2.append(x2)

# Plot them
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);
# Plot intersecting lines
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 [21]:
print(maxK1)
print(np.linalg.cond(A1))

1.0000000000000044
1.0000000000000004
In [22]:
print(maxK2)
print(np.linalg.cond(A2))

13.470236397164543
18.1448880590088
In [ ]: