HW #1

Imports and dependencies



In [1]:

    
import numpy as np

from scipy.stats import norm
from scipy.io import loadmat

Problem #1

Let $r_j$ be the $j$th row of $A$, and let $x_j$ be the solution to $$Ax_j = r_j$$



In [2]:

    
# Load the matrix into memory, assuming for now that it is stored in the home directory
A = loadmat("../dataset/hw1/A11F17108.mat")['A']

# Obtain x, where Ax[:,i] = A[i,:].T
x = np.linalg.lstsq(A,A.T)[0]



In [3]:

    
# Test
for i in range(200):
    if not np.allclose(np.dot(A,x[:,i]), A[i,:]):
        print("Fails on index %d" % i)
        break
    if i == 199:
        print("Succeeds!")









    



Succeeds!

Part (a)

Solve $$\textrm{argmin}_{j \in [199]} \|x_j\|_2$$



In [4]:

    
norms = np.linalg.norm(x,axis=0)



In [5]:

    
# Test
for i in range(200):
    if not np.allclose(np.sqrt(np.sum(np.power(x[:,i],2))), norms[i]):
        print("Fails on index %d" % i)
        break
    if i == 199:
        print("Succeeds!")









    



Succeeds!



In [6]:

    
np.argmin(norms)









    Out[6]:





164



In [7]:

    
# Test
m = np.inf
ind = 0
for i in range(200):
    if norms[i] < m:
        m = norms[i]
        ind = i
print(ind)

Part (b)

Solve $$\textrm{argmax}_{j \in [199]} \|x_j\|_2$$



In [8]:

    
np.argmax(norms)









    Out[8]:





123



In [9]:

    
# Test
m = 0
ind = 0
for i in range(200):
    if norms[i] > m:
        m = norms[i]
        ind = i
print(ind)

Part (c)

What is the average value of $\|x_j\|_2$ over all $j \in [199]$?



In [10]:

    
np.mean(norms)









    Out[10]:





47.154309366284032

Problem #2

Part (a)

Let $y$ be a vector orthogonal to $r_1^T, \dots, r_n-1^T$ with $\|y\|_2 = 1$ and $y_1 > 0$. What are $y_3$, $y_{12}$, and $y_{37}$?



In [11]:

    
# Compute the SVD of A with the last row clipped, and use the last row of V
u,s,v = np.linalg.svd(A[:-1,:])
ns = v[-1,:]



In [12]:

    
# Test
target = np.zeros((199,1))
for i in range(199):
    if not np.allclose(np.dot(A[i,:],ns),target):
        print("Fails on index %d" %i)
        break
    if i == 198:
        print("Succeeds!")









    



Succeeds!



In [13]:

    
np.linalg.norm(ns)









    Out[13]:





0.99999999999999989



In [14]:

    
# check orientation
ns[0]









    Out[14]:





0.076186565900629666



In [15]:

    
# return relevant indices
(ns[2],ns[11],ns[36])









    Out[15]:





(0.057460365244790035, 0.10102935095916442, 0.10848280840607714)



In [16]:

    
# Just in case of an oopsie
(ns[3],ns[12],ns[37])









    Out[16]:





(0.047033764877241764, 0.067764640935866607, 0.022898505380490552)



In [17]:

    
# Test
for i in range(0,199):
    if not np.allclose(np.dot(ns,A[i,:]),0):
        print("Failed at index %d" %i)
        break
    if i == 198:
        print("Succeeds!")









    



Succeeds!



In [18]:

    
q,r = np.linalg.qr(A[:-1,:].T, mode="complete")



In [19]:

    
y = q[:,-1]



In [20]:

    
y[0]









    Out[20]:





-0.07618656590063011



In [21]:

    
y = -y



In [22]:

    
(y[2],y[11],y[36])









    Out[22]:





(0.057460365244789348, 0.1010293509591651, 0.10848280840607764)



In [23]:

    
(y[3],y[12],y[37])









    Out[23]:





(0.047033764877242069, 0.067764640935865594, 0.022898505380490472)



In [24]:

    
# Test
for i in range(0,199):
    if not np.allclose(np.dot(y,A[i,:]),0):
        print("Failed at index %d" %i)
        break
    if i == 198:
        print("Succeeds!")









    



Succeeds!

Part (b)

What is the ratio $\frac{\sigma_1}{\sigma_n}$, where $\sigma_i$ is the $i$th largest singular value of $A$?



In [25]:

    
u,s,v = np.linalg.svd(A)
s[0]/s[-1]









    Out[25]:





1472.0413891117878

Problem #3

What are the first six terms of the Taylor series expansion at a point $(\alpha, \beta)$ for the function $$f(x,y) = x^2 + (y-2)^4 + x^2 y^4$$

$$\alpha^2 + (\beta-2)^2 + \alpha^2\beta^4 + 2\alpha(1+\beta^4)(x-\alpha) + 4(\alpha^2\beta^3 + (\beta-2)^3)(y-\beta)$$

Problem #4

For the same $f(x,y)$ as above, what is $\nabla f (3,1)$?

gradient

evaluation

Problem #5

Consider the curves $$C_1 = \{(x,y) \mid y = 1/x, x > 0 \}$$ $$C_2 = \{(w,z) \mid (w-5)^2 + (z-2)^2 = 1 \}$$ What points on $C_1$ and $C_2$ are closest to each other?

Consider solving $$\min_{\substack{(x,y) \in C_1 \\ (w,z) \in C_2}} d((x,y),(w,z))$$

This is equivalent to $$\min_{\substack{(x,y) \in C_1 \\ (w,z) \in C_2}} (x - w)^2 + (y - z)^2$$

Since $C_2$ defines a circle centered at $(5,2)$, a simple argument suffices to show that this is equivalent to solving $$\min_{(x,y) \in C_1} (x - 5)^2 + (y - 2)^2$$

By substitution of $y$, using the definition of $C_1$, this is equivalent to $$\min_{x > 0} \quad (x - 5)^2 + (x^{-1} - 2 )^2$$

The derivative of this function is $$2(-x^{-3} + 2x^{-2} + x -5)$$

Solve for x when set equal to zero



In [26]:

    
# solve for y
x = 4.92594
y = 1./x



In [27]:

    
y









    Out[27]:





0.2030069387771674



In [28]:

    
(x,y)









    Out[28]:





(4.92594, 0.2030069387771674)



In [29]:

    
# solve for the vector connecting (x,y) to (5,2)
s = np.array([x-5,y-2])
s









    Out[29]:





array([-0.07406   , -1.79699306])



In [30]:

    
# normalize and scale
nm = np.linalg.norm(s)
w = s[0]/nm+5
z = s[1]/nm+2
(w,z)









    Out[30]:





(4.9588216644570524, 1.0008481873699511)



In [31]:

    
# Test
(w-5)**2 + (z-2)**2









    Out[31]:





1.0



In [32]:

    
# Print the distance
np.linalg.norm([x-w,y-z])









    Out[32]:





0.79851854193472449

Of course, this solution can also be obtained by way of stochastic gradient descent (not implemented here).

Problem #6

Shuffle a normal deck of 52 playing cards and draw three cards.

Part (a)

What is the probability that they are the same suit in increasing rank?



In [33]:

    
def choose(n, k):
    """
    A fast way to calculate binomial coefficients by Andrew Dalke (contrib).
    """
    if 0 <= k <= n:
        ntok = 1
        ktok = 1
        for t in xrange(1, min(k, n - k) + 1):
            ntok *= n
            ktok *= t
            n -= 1
        return ntok // ktok
    else:
        return 0



In [34]:

    
# the total number of possible triples
tot = choose(52,3)*6



In [35]:

    
# the total number of acceptable triples over the total number of triples
float(choose(13,3)*4)/tot









    Out[35]:





0.008627450980392156



In [36]:

    
# Test
# Count all possible acceptable triples and divide by all possible triples
trips = [(i,j,k) for i in range(1,12) for j in range(i+1,13) for k in range(j+1,14)]
float(len(trips)*4)/tot









    Out[36]:





0.008627450980392156



In [37]:

    
# Examine the acceptable triples. Bear in mind there are four of each type
trips[:20]









    Out[37]:





[(1, 2, 3),
 (1, 2, 4),
 (1, 2, 5),
 (1, 2, 6),
 (1, 2, 7),
 (1, 2, 8),
 (1, 2, 9),
 (1, 2, 10),
 (1, 2, 11),
 (1, 2, 12),
 (1, 2, 13),
 (1, 3, 4),
 (1, 3, 5),
 (1, 3, 6),
 (1, 3, 7),
 (1, 3, 8),
 (1, 3, 9),
 (1, 3, 10),
 (1, 3, 11),
 (1, 3, 12)]

Part (b)

Given that the first two cards you drew from the deck are of different suits, what is the probability that they are all of different suits?



In [38]:

    
# Much simpler.
26./50









    Out[38]:





0.52

Problem #7

The random variables $X$ and $Y$ have a joint pdf $p(X,Y)$ with distribution $(X,Y) \sim \mathcal{N}((\mu_X,\mu_Y)),\Sigma)$, where $(\mu_X,\mu_Y) = (2, -1)$ and $$ \Sigma = \begin{bmatrix} 1 & 0 \\ 0 & 2 \end{bmatrix}$$

Assume throughout that $(x,y) \sim p(X,Y)$

Part (a)

What is $\mathbb{P} ( x \in [0,2])$?

Note that the covariance matrix of $p(X,Y)$ has zeros on the cross-diagonal entries. This implies that $X$ and $Y$ are independent of one another. Hence, we have that $X \sim \mathcal{N}(2,1)$ and $Y \sim \mathcal{N}(-1,2)$.

Thus, $\mathbb{P} ( x \in [0,2]) = \int_0^2 \int_{-\infty}^\infty p(X,Y) dY dX = \int_0^2 p(X) dX$



In [39]:

    
# Create given normal variables
X = norm(loc=2,scale=1)
Y = norm(loc=-1,scale=np.sqrt(2))



In [40]:

    
X.cdf(2) - X.cdf(0)









    Out[40]:





0.47724986805182079



In [41]:

    
Y.cdf(2) - Y.cdf(0)









    Out[41]:





0.22280263433113212

Part (b)

What is $\mathbb{P} ( x \in [0,2] \vee y \in [0,2])$?



In [42]:

    
(X.cdf(2) - X.cdf(0)) + (Y.cdf(2) - Y.cdf(0)) - ((X.cdf(2) - X.cdf(0))*(Y.cdf(2) - Y.cdf(0)))









    Out[42]:





0.59371997454682202

Part (c)

What is $\mathbb{P} ( x \in [0,2] \wedge y \in [0,2])$?



In [43]:

    
(X.cdf(2) - X.cdf(0))*(Y.cdf(2) - Y.cdf(0))









    Out[43]:





0.10633252783613088