Python Numpy.array() Value error as a rite of passage?

[MikeLawrence]

Is it a rite of passage to type up a NLA algo using numpy and then

get hit with a dozen errors stimming from incorrect shape configuration in matrix multiplication​

search for an hour or two wondering how they came to be and snuffing them all out​

use np.transpose in many random places​

finally realize that np.array() defaults to a 1xn vector but standard math notation uses nx1 vectors?​

​

Because that's what I just did, and I'd like to think I'm not alone.

As a follow up, is it standard to transpose all vectors upon input or creation or does everyone just accept the flipped dimensions and move on. I assume we accept and move on.

 

[MikeLawrence]

I've been using np.array's for a bit now, and I've done a tiny bit of NLA in C++ for the final project of the first course, but I hadn't coded in a bit and for some reason I either just forgot that the default was 1xn or it wasn't an issue in whatever I've been doing.

All the notation, and also how I was thinking about vectors, was that it was just 'n dimensional.' I suppose I assumed that python/numpy would handle it smoothly depending on what you were multiplying it against.

The 'aha!' moment was very satisfying, but man do I feel slow.

In a related note, we can multiply constants into vectors, but only on the side that has one dimension. Why is that? It's a constant. Python knows it's a constant. The error message I get tells me it is 1x1. BUT, it is multiplied like a vector.
edit: I actually think this one is due to my 'constant' being the result of the (1x10)(10x1) dot product, so it actually is held as a 1x1 vector. I'll leave my confusion in for future confused python learners.
edit#2: yeah that's got to be why, I double checked and you definitely can multiply constants on any side.

 

[Paul Lopez]

rite of passage absolutely

When you multiply a constant by a vector (or a matrix) in Python, you typically use a broadcasting mechanism. Broadcasting is a method that allows you to perform operations on arrays of different shapes by automatically expanding the smaller array to match the shape of the larger one. Def read up on that when you can. It'll help with these lovely yet annoying things that happen in python.

 

[Daniel Duffy]

Take a step backwards and focus on the problem (ignore the bells and whistles).
numpy is only a library.

#QR decomposition
r = c = 4
A = np.random.randn(r, c); print(A)
Q,R = np.linalg.qr(A); print(Q); print(R)
# Cross-check
print (np.matmul(Q,R))

QR = np.matmul(Q,R)
for i in range(0,r):
    for j in range(0,c):
        A[i,j] = A[i,j] - QR[i,j]
print(A) # should be very much zero

#svd
A = np.random.randn(9,6) + 1j* np.random.randn(9,6)
U,s,V = np.linalg.svd(A, full_matrices = True)
print(U.shape, V.shape, s.shape) #(9,9), (6,6), (6,)

S = np.zeros((9,6), dtype = complex) #S = [[s,0],[0.0]]
S[:6,:6]=  np.diag(s) #embed vector in matrix

pi = np.dot(U, np.dot(S,V)) #pseudoinverse
print(np.allclose(A, pi)) #function to compare two arrays

 

[Daniel Duffy]

finally realize that np.array() defaults to a 1xn vector but standard math notation uses nx1 vectors?

Maybe avoid defaults; ReadTheDocs?

import scipy.sparse.linalg as spla
import numpy as np
from numpy import linalg as LA
from scipy import linalg
import random

print("Solve Ax = b")
A = np.array([[3, 1],[1,2]]); b = np.array([9,8])
x = np.linalg.solve(A,b); print(x) # (2,3)

#Conjugate gradient DD 2019-12-14
y = spla.cg(A, b); print(y) # (2,3)

# Wiki example
A = np.array([[4, 1],[1,2]]); b = np.array([1,2])
z = np.linalg.solve(A,b); print(z) # (0.0909, 0.6364)

 

[Daniel Duffy]

rite of passage absolutely

When you multiply a constant by a vector (or a matrix) in Python, you typically use a broadcasting mechanism. Broadcasting is a method that allows you to perform operations on arrays of different shapes by automatically expanding the smaller array to match the shape of the larger one. Def read up on that when you can. It'll help with these lovely yet annoying things that happen in python.

Rule of thumb

Let M(n,m) == matrices with n rows and m columns

A, B matrices then C = A*B is possible if #columns(A) = #rows(B)

A ~ M(n,p), B ~= M(p,m)
C = A*B

then C ~ M(n,m).

also applies when B is a vector.

 

[MikeLawrence]

finally realize that np.array() defaults to a 1xn vector but standard math notation uses nx1 vectors?

Maybe avoid defaults; ReadTheDocs?

import scipy.sparse.linalg as spla
import numpy as np
from numpy import linalg as LA
from scipy import linalg
import random

print("Solve Ax = b")
A = np.array([[3, 1],[1,2]]); b = np.array([9,8])
x = np.linalg.solve(A,b); print(x) # (2,3)

#Conjugate gradient DD 2019-12-14
y = spla.cg(A, b); print(y) # (2,3)

# Wiki example
A = np.array([[4, 1],[1,2]]); b = np.array([1,2])
z = np.linalg.solve(A,b); print(z) # (0.0909, 0.6364)

print(A.shape, x.shape, b.shape)
#Output: (2,2), (2,),(2,)

I've done it like this as well, and it does work, but it's weird in that python doesn't know what the second dimension of the arrays is.
I thought that left multiplying with such a (n,?) vector gives an error, but I'm testing it again now and it appears this is the perfect workaround.

I'm loosing my mind, I can't get np.arrays to default to (1,n) vectors anymore. (array name).shape is always returning (n,). What was I doing before??

Thanks anyway Dr. Duffy.

 

[Quasar Chunawala]

If you pass in a list in R^n to numpy.array(), it always creates an object of shape (n,). If you pass in m lists, each consisting of n lists of p elements, i.e. something in R^(m×n×p) to numpy.array(), it always creates an object of shape (m,n,p). This stuff is pretty well documented.

A numpy.array of shape (1,1) is still a matrix at the end of the day(and not a python int/float). You could np.matmul() or np.dot() it with another object as long as you follow the usual matrix multiplication rules. (For np.dot(A,B), numpy internally checks if the condition A.shape[-1] == B.shape[0] is met)

There is no obligation to transpose a numpy array. It depends on the use-case. If A.shape = (3,3) and x.shape=(3,), np.dot(A,x) should work out of the box.

 

[Quasar Chunawala]

One fun exercise - to get your hands wet with numpy vectorization. Generate paths of your favorite stochastic process such as OU/Brownian Bridge.

Step 1. Draw standard gaussians of size (num_steps,num_paths). Let's call this Z. These are Brownian increments.
Step 2. Construct the covariance matrix C of your process.
Step 3. Perform Cholesky factorization of C as C=(A A^T)
Step 4. Multiply X=AZ to get increments of your process. X will have shape (num_steps, num_paths).
Step 5. Add up these increments to get all paths of your favorite process.

Obviously, there should be no loops (except of step 2)!