Question a#

Carry out the matrix-vector multiplication $Axx$ in SymPy, NumPy and Maple. Hint: The matrix-vector and thus matrix-matrix multiplication is done using different syntax in all three programs, which will be using either * (asterisk), . (full-stop/period) or @ (at).

Question b#

The expression $A+2$ has no mathematical definition – so it makes no mathematical sense. Investigate whether it makes sense in SymPy, NumPy, and Maple, or whether you are given an error.

Again all three programs SymPy, NumPy, and Maple behave differently. Discuss which program you think behaves the most reasonable.

Question c#

NumPy often behaves like SymPy (and Maple), for instance when transposing a matrix:

A.T, A_np.T

(Matrix([
 [1, 4, 7],
 [2, 5, 8],
 [3, 6, 9]]),
 array([[1., 4., 7.],
        [2., 5., 8.],
        [3., 6., 9.]]))

but NumPy can create many other smart manipulations with arrays. This also means that one has to be careful not to make mistakes. It is important to know that NumPy does not know the difference between row- and column-vectors: All vectors in NumPy are 1-dimensional. One can thus not say whether it is a “upright” vector or a “lying-down” vector, so this distinction would make no difference in NumPy. This can be seen by x_np having a size of “3” rather than “(3,1)” or “(1,3)”:

A_np.shape, x_np.shape

In the following you must explain how/what NumPy computes, and how you would write out mathematically or in SymPy (if at all possible):

x_np @ x_np, x_np.dot(x_np)

A_np.dot(x_np)

array([14., 32., 50.])

x_np @ A_np

array([30., 36., 42.])

A_np * x_np, x_np * A_np

(array([[ 1.,  4.,  9.],
        [ 4., 10., 18.],
        [ 7., 16., 27.]]),
 array([[ 1.,  4.,  9.],
        [ 4., 10., 18.],
        [ 7., 16., 27.]]))

Question d#

When we create orthogonal projections on lines as in Week 3, we often need to create an $n\times n$ projection matrix from the matrix-matrix product of an $n\times 1$ column vector and a $1\times n$ row vectors. In SymPy we would do:

u = x / x.norm()
u * u.T

_.evalf()   # as a decimal number

How can you do this in NumPy? Why does the following not work?

u_np = x_np / np.linalg.norm(x_np)
u_np @ u_np.T

Question e#

Download the file weight-height.csv from Learn. The data file is from Kaggle and can also be downloaded from there. Import the CSV file with NumPy using the following command

# Load the CSV file into a NumPy array
data = np.genfromtxt('weight-height.csv', delimiter=',', dtype=[('Gender', 'U10'), ('Height', float), ('Weight', float)], names=True)

# Access the columns by their names
gender = data['Gender']
height = data['Height']
weight = data['Weight']

# Print the data
# print(data.dtype)
print(gender)
print(height)
print(weight)

['"Male"' '"Male"' '"Male"' ... '"Female"' '"Female"' '"Female"']
[73.84701702 68.78190405 74.11010539 ... 63.86799221 69.03424313
 61.94424588]
[241.89356318 162.31047252 212.74085556 ... 128.47531878 163.85246135
 113.64910268]

Question f#

Explain the following plot. What does each point/intersection correspond to?

plt.scatter(x=height, y=weight, color='c', marker='x')
plt.title("Height versus Weight")
plt.xlabel("Height (inches)")
plt.ylabel("Weight (lbs)")
plt.plot() 

Question g#

Explain what the following data matrix contains:

X = np.array([height, weight]).T
X.shape

# the first 10 rows
X[:,0:10] # or just X[:10]

array([[ 73.84701702, 241.89356318],
       [ 68.78190405, 162.31047252],
       [ 74.11010539, 212.74085556],
       ...,
       [ 63.86799221, 128.47531878],
       [ 69.03424313, 163.85246135],
       [ 61.94424588, 113.64910268]])

Question h#

Write out a Python function that computes the average of each column in the $X$ matrix:

def average_of_each_column(X):
    # add code
    return X[0,:]  # fix this to return the average of each column

The output (after return) must be a NumPy array with shape(2,), meaning an average-value vector with a length fo 2. Check your function with the call:

X - average_of_each_columns(X)

Question i#

We must now standardize the data, which in this exercise is done simply by centering the data set (the “data cloud”) at $(0,0)$. More precisely, in each column of $X$ you must subtract the column average from each element in the column. In Python simply:

def standardize(X):
    return X - average_of_each_column(X) # 'broadcasting' the vector as we saw above

X_st = standardize(X)

Plot the new standardized data set $X_{st}$ and check that it is centered around $(0,0)$.

Hint: Remember that the data set $X$ can be plotted by:

plt.scatter(X[:,0], X[:,1]) # X[:,0] is height and X[:,1] is weight
plt.title("Height versus Weight")
plt.xlabel("Height (inches)")
plt.ylabel("Weight (lbs)")
plt.plot() 

Now do the same for $X_{st}$.

Question j#

We will now briefly return to the mathematics for a moment and ask you to prove the following:

Let $X$ be an arbitrary, real $n\times k$ matrix. Proof that the $k\times k$ matrix $C$ given by $C=X^{T}X$ is symmetric. Argue that $C$ has real eigenvalues and $k$ orthogonal eigenvectors.

Question k (optional)#

Prove that $C$ actually is positive semi-definite and thus cannot have negative eigenvalues. You may skip this exercise for now and return to it later.

Question l#

We return to the standardized data matrix $X_{st}$ and compute the $C$ matrix by:

C = 1/(X_st.shape[0]-1) * X_st.T @ X_st
C

array([[  14.80347264,  114.24265645],
       [ 114.24265645, 1030.95185544]])

The matrix $C$ is called the (sample) covariance matrix. The constant 1/(X_st.shape[0]-1), which is $1/(n-1)$, is not important to us. But we are interested in the eigenvalues and eigenvectors of $C$ which are found with the call:

lamda, Q = np.linalg.eig(C)
lamda, Q

(array([   2.11786479, 1043.63746329]),
 array([[-0.99389139, -0.1103626 ],
        [ 0.1103626 , -0.99389139]]))

Show that $Q$ already has normalized columns. From this it follows that $Q$ is (real) orthogonal. Why?

Question m#

Create a $2\times 2$ projection matrix $P$ that describes the orthogonal projection down onto the the eigenvector space corresponding to the largest eigenvalue of $C$. Let $x{x}_k$ be the $k$’th row in $X_{st}$, meaning the weight and height of the $k$’th test subject. As a column vector this is written as $(x{x}_k{)}^T$, which we can project down onto the eigenvector space corresponding to the largest eigenvalue by $P(x{x}_k{)}^T$. If we wish to do this for all test subjects, we thus just have to compute $P(X_{st}{)}^T$. The matrix $X_{st}$ is of size $n\times 2$, while $P(X_{st}{)}^T$ is of $2\times n$, so we transpose this in order to achieve the same size as $X_{st}$, all together $(P(X_{st}{)}^T{)}^T=X_{st}{P}^T$.

First, create the $P$ matrix. Remember that you can extract e.g. the second column of $Q$ with q2 = Q[:,1]. Check that $P$ is a $2\times 2$ projection matrix. Then plot the standardized data set $X_{st}$ along with the projected data set $X_{st}{P}^T$. You must end with a plot that looks as follows:

About Python plots: You can plot subspaces spanned by columns of a matrix with the following code:

q1 = Q[:,0] # the first column
q2 = Q[:,1] # the second column

# Define a range for t, which will be used to extend the lines
t = np.linspace(-10, 10, 2)

# Plotting
fig, ax = plt.subplots()
# For each vector, plot a line that it spans
ax.plot(t*q1[0], t*q1[1], 'r', label='Line Spanned by Column 1')
ax.plot(t*q2[0], t*q2[1], 'b', label='Line Spanned by Column 2')

# Adjust the plot limits and style
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.grid(color = 'gray', linestyle = '--', linewidth = 0.5)
plt.legend()
plt.title('Lines Spanned by the Columns of Q')

# Display the plot
plt.show()

Our data set $X$ (and $X_{st}$) has two components, those being 2 eigenvalues/eigenvectors. Via the projection $X_{st}{P}^T$ we carry out a dimensional reduction (from 2 to 1) down onto the most important dimension of the data set. In this example it is not smart to do this dimensional reduction, since the variance in both weight and height are important, but we have made the best representation if we only have one “component/coordinate vector”.

If the data $X_{st}$ had $k$ columns and the $C$ matrix thus became $k\times k$ with $k$ eigenvalues, one should do an orthogonal projection of $X_{st}$ onto the subspace of ${\mathbb{R}}^k$ spanned by the eigenvectors corresponding to the “largest” eigenvalues, since these are the “most important” components. If the eigenvalues of $C$ for instance are $1,0.7,0.5,0.01,001$, $P$ should do the orthogonal projection on the subspace spanned by the eigenvectors corresponding to the three largest eigenvalues.

Question n#

Let us again look at some data matrices $Y$ with 3 columns such that we can consider $Y$ as a set of points in ${\mathbb{R}}^3$. The data matrix $Y$ in the code below is generated “randomly”, but according to the rules:

1. $Y$ is standardized (each column has zero as its average),

2. The eigenvalues of $C=\frac{1}{n-1}{Y}^{T}Y$ are approximately equal to ${\lambda }_1,{\lambda }_2.{\lambda }_3$, which are given by the use in the line eigenvalues = [5, 1, 0.04].

# Function to generate a synthetic dataset
def generate_3d_data(eigenvalues, size=1000, dim=3):
    # Eigenvalues specified by the user
    assert len(eigenvalues) == dim, "There must be exactly dim eigenvalues."
    
    # Create a diagonal matrix for the eigenvalues
    Lambda = np.diag(eigenvalues)
    
    # Generate a random orthogonal matrix (eigenvectors)
    Q, _ = np.linalg.qr(np.random.randn(dim, dim))
    
    # Generate random data from a standard normal distribution
    data_standard_normal = np.random.randn(size, dim)
    
    # Transform the data using the square root of the covariance matrix
    data_transformed = data_standard_normal @ np.sqrt(Lambda) @ Q.T
        
    # Construct the covariance matrix
    Cov = 1/(size-1) * data_transformed.T @ data_transformed
     
    return data_transformed, Cov

# Eigenvalues you want for your covariance matrix
eigenvalues = [5, 1, 0.04]

# Generate the data
Y, C = generate_3d_data(eigenvalues)

We can plot the data matrix $Y\in {\mathbb{R}}^{1000\times 3}$ as 1000 data points in ${\mathbb{R}}^3$:

# %matplotlib qt   # plotted in a new window which can be rotated (remove #)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(Y[:,0], Y[:,1], Y[:,2], c='b', marker='x')
ax.set_xlabel('first column') 
ax.set_ylabel('second column')
ax.set_zlabel('third column')
plt.title('3D Scatter Plot of Synthetic Data')
plt.show()

Choose different positive ${\lambda }_1,{\lambda }_2,{\lambda }_3$ in eigenvalues = [5, 1, 0.04] and investigate the appearing set of points by rotating the plot (you must remove the comment symbol from %matplotlib qt). Try both with examples where ${\lambda }_1$ is much larger than ${\lambda }_2,{\lambda }_3$, and where ${\lambda }_3$ is much smaller than ${\lambda }_1$ and ${\lambda }_2$. How does the shape of the set of points change with the eigenvalues?

Think about this (roughly formulated statement): Eigenvalues of $C$ are, roughly speaking, given by the semi-axes in the ellipsoid that most effectively can “cover” the most data points. The directions of the symmetry axes of the ellipsoid are determined by the eigenvectors. Non-important “components” are along direction in which the ellipsoid is “flat”.

Question o#

In this Theme exercise we are centering the dataset by subtracting the mean value from each column. It can also be useful (or necessary) to “scale” the data, such that all columns get the same spread (meaning, the variation in each column becomes the same). This standardization is called z-score normalization, where the data are transformed such that each column get the mean value 0 and spread 1.

Discuss when it can be necessary to standardize the data in this way (before the covariance matrix is calculated and PCA is performed).

Question p#

Perform PCA on one or more of the synthetic data sets $Y$ in ${\mathbb{R}}^3$ from the previous problem.

Question q#

Return to the data set $X$. Perform PCA for, respectively, male and female test subjects:

# Filter the data for males
male_data = data[data['Gender'] == '"Male"']
female_data = data[data['Gender'] == '"Female"']

# Access the columns for male data
male_height = male_data['Height']
male_weight = male_data['Weight']
female_height = female_data['Height']
female_weight = female_data['Weight']

# Plot the male weight and height
plt.scatter(x=male_height, y=male_weight, color='c', label='Male')
plt.scatter(x=female_height, y=female_weight, color='b', label='Female')
plt.title("Female vs Male Height and Weight")
plt.xlabel("Height (inches)")
plt.ylabel("Weight (lbs)")
plt.legend()
plt.plot()