Stanford CS205a¶

Norms, Sensitivity & Conditioning¶

Supplemental Julia Notebook

Prof. Doug James

Vector p-norms in Julia¶

In [1]:

?(norm)

search: norm normpath normalize normalize! normalize_string vecnorm issubnormal

Out[1]:

norm(A, [p])

Compute the p-norm of a vector or the operator norm of a matrix A, defaulting to the p=2-norm.

For vectors, p can assume any numeric value (even though not all values produce a mathematically valid vector norm). In particular, norm(A, Inf) returns the largest value in abs(A), whereas norm(A, -Inf) returns the smallest.

For matrices, the matrix norm induced by the vector p-norm is used, where valid values of p are 1, 2, or Inf. (Note that for sparse matrices, p=2 is currently not implemented.) Use vecnorm to compute the Frobenius norm.

In [2]:

a=ones(3)

Out[2]:

3-element Array{Float64,1}:
 1.0
 1.0
 1.0

In [3]:

norm(a)

Out[3]:

1.7320508075688772

In [4]:

norm(a,2)

Out[4]:

1.7320508075688772

In [5]:

norm(a,Inf)

Out[5]:

1.0

In [6]:

norm(a,1)

Out[6]:

3.0

Visualizing Vector p-norms¶

Let us look at contour plots of \begin{equation} \| (x,y) \|_p \end{equation} for different values of $p$.

In [7]:

using Plots

In [8]:

f(x,y) = vecnorm([x y], p)   # parameter-dependent p-norm function

Out[8]:

f (generic function with 1 method)

In [9]:

x = y = linspace(-1, 1, 100);

In [10]:

p=1;  
contourf(x, y, f, xlabel="x", ylabel="y", title=string("Norm, p=",p))

Out[10]:

In [11]:

p=2;  
contourf(x, y, f, xlabel="x", ylabel="y", title=string("Norm, p=",p))

Out[11]:

In [12]:

p=5;  
contourf(x, y, f, xlabel="x", ylabel="y", title=string("Norm, p=",p))

Out[12]:

In [13]:

p=Inf;  
contourf(x, y, f, xlabel="x", ylabel="y", title=string("Norm, p=",p))

Out[13]:

Matrix-transforming vectors on a circle¶

In [15]:

unitVec(θ)=[cos(θ), sin(θ)]  # Unit vector (in 2-norm) in direction θ

WARNING: Method definition unitVec(Any) in module Main at In[14]:1 overwritten at In[15]:1.

Out[15]:

unitVec (generic function with 1 method)

In [16]:

unitVec(0)

Out[16]:

2-element Array{Float64,1}:
 1.0
 0.0

In [17]:

UV=zeros(2,360)   # Matrix of unit 2-vectors as columns
for i=1:360
    UV[:,i] = unitVec(i/π)
end
UV

Out[17]:

2×360 Array{Float64,2}:
 0.949766  0.80411   0.577666  0.293185  …  0.654252  0.3847    0.076498
 0.312962  0.594481  0.816273  0.956056     0.756277  0.923042  0.99707

In [29]:

R = randn(2,2) # randn matrix
B = R * UV     # Multiply UV by random matrix
plot(B[1,:], B[2,:], title=string("RandnMtx * UnitCircle, ||R||_2=",norm(R,2)))

Out[29]:

2D Visualization of Induced Matrix Norms¶

For a 2D unit vector $\vec{x}(\theta) = (\cos\theta, \sin\theta)$, and a random matrix $\mathbf{R}$, let's visualize the induced matrix $p$-norm ratios \begin{equation} ratio(\theta)=\frac{\| \mathbf{R} \; \vec{x}(\theta) \|_p}{\|\vec{x}(\theta)\|_p} \end{equation} for different values of $p$.

In [30]:

nRUVRatio(θ) = norm(R*unitVec(θ), p)/norm(unitVec(θ),p)  # p-norm of R*unitVec / p-norm of unitVec

Out[30]:

nRUVRatio (generic function with 1 method)

In [31]:

R=randn(2,2)  # a randn matrix to measure
println("R = ",R)

R = [0.504342 0.139108; -0.0371534 -0.0523837]

In [32]:

p=Inf;
println("||R||_",p," = ",norm(R,p))
plot(nRUVRatio, linspace(0,π,360), xlabel="θ", title=string("ratio(θ) for p=",p))

||R||_Inf = 0.6434504854188444

Out[32]:

In [33]:

p=2;
println("||R||_",p," = ",norm(R,p))
plot(nRUVRatio, linspace(0,π,360), xlabel="θ", title=string("ratio(θ) for p=",p))

||R||_2 = 0.5255487353898564

Out[33]:

In [34]:

p=1;
println("||R||_",p," = ",norm(R,p))
plot(nRUVRatio, linspace(0,π,360), xlabel="θ", title=string("ratio(θ) for p=",p))

||R||_1 = 0.5414954937506968

Out[34]:

Experiments with an ill-conditioned Vandermonde matrix¶

In [35]:

# A famous ill-conditioned matrix: Vandermonde matrix
vander(N) = [x^i for x=linspace(1./N,1,N), i=1:N] # Vandermonde matrix

Out[35]:

vander (generic function with 1 method)

In [36]:

N = 4
V = vander(N)

Out[36]:

4×4 Array{Float64,2}:
 0.25  0.0625  0.015625  0.00390625
 0.5   0.25    0.125     0.0625    
 0.75  0.5625  0.421875  0.316406  
 1.0   1.0     1.0       1.0

In [37]:

println("size = ",size(V))
plot(V, title=string("Vandermonde matrix. size=",size(V)))

size = (4,4)

Out[37]:

In [38]:

cond(V)   # Condition number (expensive computation but...)

Out[38]:

524.7924487965224

In [39]:

svdvals(V)

Out[39]:

4-element Array{Float64,1}:
 2.3134    
 0.447495  
 0.0601855 
 0.00440821

In [42]:

b=ones(N)

Out[42]:

4-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0

In [43]:

LU = lufact(V)
x=LU\b

Out[43]:

4-element Array{Float64,1}:
   8.33333
 -23.3333 
  26.6667 
 -10.6667

In [44]:

vecnorm(b-V*x)/vecnorm(b)   # relative error of matrix solve

Out[44]:

8.899114524108741e-16

How does solve rel-error depend on size of Vandermonde matrix?¶

In [45]:

function vanderSolveRelErr(N)
    V = [x^i for x=linspace(1.0/N,1.0,N), i=1:N]
    b = ones(N)  # a predictable RHS, \vec{1}
    LU = lufact(V)
    x = LU \ b
    return vecnorm(b-V*x)/vecnorm(b)
end

Out[45]:

vanderSolveRelErr (generic function with 1 method)

In [46]:

vErr = [vanderSolveRelErr(n) for n=1:100]   # solve all problems from size n=1 to 100

Out[46]:

100-element Array{Float64,1}:
 0.0        
 0.0        
 2.56395e-16
 8.89911e-16
 4.19071e-15
 1.07055e-14
 3.29288e-14
 1.09224e-13
 3.2904e-13 
 2.90792e-12
 1.02942e-11
 4.54925e-11
 1.48135e-10
 ⋮          
 0.0906789  
 1.94849    
 0.198944   
 0.363464   
 3.10338    
 0.321879   
 0.299673   
 0.614996   
 0.0526596  
 0.0704173  
 0.0933069  
 0.306217

In [47]:

plot(log10.(vErr+eps()), title="RelError solving V x = 1")

Out[47]:

Testing $V^{-1} V = I$:

In [56]:

for n=2:2:50
    println("n=",n,": relErr(invV*V-I): ",norm(inv(vander(n))*vander(n)-eye(n))/norm(eye(n)))
end

n=2: relErr(invV*V-I): 0.0
n=4: relErr(invV*V-I): 5.704378753739517e-15
n=6: relErr(invV*V-I): 7.206346992435156e-13
n=8: relErr(invV*V-I): 4.232127554641132e-11
n=10: relErr(invV*V-I): 3.494164228312279e-9
n=12: relErr(invV*V-I): 2.67328213816742e-7
n=14: relErr(invV*V-I): 1.13656866037886e-5
n=16: relErr(invV*V-I): 0.0005713609123286228
n=18: relErr(invV*V-I): 0.03460692390723406
n=20: relErr(invV*V-I): 2.7953035818633856
n=22: relErr(invV*V-I): 32.06809436167812
n=24: relErr(invV*V-I): 2981.979751536043
n=26: relErr(invV*V-I): 74.36798538478568
n=28: relErr(invV*V-I): 146.86182945261046
n=30: relErr(invV*V-I): 61.98487651414396
n=32: relErr(invV*V-I): 632.6816419885412
n=34: relErr(invV*V-I): 980.6993965209957
n=36: relErr(invV*V-I): 86.15079649282958
n=38: relErr(invV*V-I): 205.17510463084815
n=40: relErr(invV*V-I): 1207.7221166111324
n=42: relErr(invV*V-I): 215.92511095558865
n=44: relErr(invV*V-I): 1380.6409470134085
n=46: relErr(invV*V-I): 170.36879116352844
n=48: relErr(invV*V-I): 867.0948304744451
n=50: relErr(invV*V-I): 286.17763388790155

PROBLEM: Vandermonde matrix has many near-zero singular values for large $N$.

In [51]:

V = vander(100)
plot(log10.(svdvals(V)))

Out[51]:

In [52]:

det(V)   # numerically singular

Out[52]:

-0.0

In [53]:

cond(V)  # awful condition number

Out[53]:

5.205505037741025e19

In [ ]: