BDS 761: Data Science and Machine Learning I
Topic 4: Precision and Norms
This topic:¶
- Numerical precision (on modern desktop computers)
- Norms, distances, and basic stats
- k-means clustering
Readings:
Numerical Precision¶
Binary numbers¶
https://mikkegoes.com/computer-science-binary-code-explained/
The basis for all modern computing.
What set of numbers can be represented with the above binary scheme and 8 bits?
What is the largest number possible with 8 bits? with 32 bits?
How many decimal digits do these amount to roughly?
bin(216) # gives binary representation of integer
'0b11011000'
import numpy as np
np.base_repr(216, base=2, padding=0)
'11011000'
Hexadecimal¶
Base-16 system. Effectively groups four binary digits into each hex digit.
for k in range(0,32):
print(f'{k:4d} {bin(k):>8s} {hex(k):>5s}.')
0 0b0 0x0. 1 0b1 0x1. 2 0b10 0x2. 3 0b11 0x3. 4 0b100 0x4. 5 0b101 0x5. 6 0b110 0x6. 7 0b111 0x7. 8 0b1000 0x8. 9 0b1001 0x9. 10 0b1010 0xa. 11 0b1011 0xb. 12 0b1100 0xc. 13 0b1101 0xd. 14 0b1110 0xe. 15 0b1111 0xf. 16 0b10000 0x10. 17 0b10001 0x11. 18 0b10010 0x12. 19 0b10011 0x13. 20 0b10100 0x14. 21 0b10101 0x15. 22 0b10110 0x16. 23 0b10111 0x17. 24 0b11000 0x18. 25 0b11001 0x19. 26 0b11010 0x1a. 27 0b11011 0x1b. 28 0b11100 0x1c. 29 0b11101 0x1d. 30 0b11110 0x1e. 31 0b11111 0x1f.
Floating point numbers¶
IEEE 754 standard for representing a real number with a 32-bit (or 64-bit) word.
$$x = (-1)^{\text{sign}} \times \text{fraction} \times 2^\text{exponent}$$
Similar approach to achieve higher precision using 64-bit words
Note that two version of zero are possible, +0 and -0, depending on sign bit
$\pm$ Infinity = exponents all 1's, fraction all 0's
Not-a-number (NaN) = exponents all 1's, fraction not all 0's.
https://learn.microsoft.com/en-us/cpp/build/ieee-floating-point-representation?view=msvc-170
1./0.
--------------------------------------------------------------------------- ZeroDivisionError Traceback (most recent call last) ~\AppData\Local\Temp\ipykernel_26264\832528835.py in <module> ----> 1 1./0. ZeroDivisionError: float division by zero
import numpy as np
1/np.array([0.]), -1./np.array([0.])
C:\Users\micro\AppData\Local\Temp\ipykernel_26264\3923875365.py:3: RuntimeWarning: divide by zero encountered in true_divide 1/np.array([0.]), -1./np.array([0.])
(array([inf]), array([-inf]))
0.0/0.0
--------------------------------------------------------------------------- ZeroDivisionError Traceback (most recent call last) ~\AppData\Local\Temp\ipykernel_26264\1949058463.py in <module> ----> 1 0.0/0.0 ZeroDivisionError: float division by zero
0.0/np.array([0.0])
C:\Users\micro\AppData\Local\Temp\ipykernel_26264\844779733.py:1: RuntimeWarning: invalid value encountered in true_divide 0.0/np.array([0.0])
array([nan])
0*np.nan
nan
np.inf*np.nan
nan
Commonly used number types¶
| Common names | bits |
|---|---|
| character (char) | 8 |
| integers (int) | 32 |
| unsigned integers (uint) | 32 |
| floating-pointing numbers (float, single) | 32,64 |
| double precision floating point numbers (double) | 64 |
| short integers (short) | 16 |
| long integers (long) | 32 |
| double long integers (longlong) | 64 |
| unsigned long integers (ulong) | 64 |
Note that number of bits used vary depending on hardware, compiler, ...
Precision (a.k.a. number of bits) is often added to the type name, e.g. float32.
Working with limited precision¶
Number of decimal digits from the number of binary exponent bits
- 32 bit floats have about 7 digits of precision
- 64 bit floats have about 16 digits precision
Must keep in mind that decimal numbers are always approximate.
1 + 1e-12
1.000000000001
1 + 1e-25
1.0
1 + 1e-25 == 1
True
1e-25 + 1e-12
1.0000000000001e-12
A = np.random.randint(0, 10, size=(3, 3))
A_inv = np.linalg.inv(A)
print(A@A_inv) # matrix times its inverse should be identity
[[ 1.00000000e+00 0.00000000e+00 -2.22044605e-16] [-3.70074342e-18 1.00000000e+00 5.55111512e-17] [ 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
Floating-Point Arithmetic¶
Floating point computations are not exact
$$ a+b \approx c$$Result will be nearest floating point number that can be represented with the binary system used.
hard zeros ~ special case when "0.0" really is all zeros (for fractional part) in binary
Logic, however, is still exact
0.1 + 0.1 + 0.1 == 0.3
False
0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1.0
False
round(0.1 + 0.1 + 0.1, 5) == round(0.3, 5) # round to 5 digits precision
True
import math
math.isclose(0.1 + 0.1 + 0.1, 0.3)
True
Overflow and Underflow¶
import numpy as np
print("\n--- Integer Overflow (int8) ---")
a = np.int8(120)
b = np.int8(10)
overflow_result = a + b
print("a =", a)
print("b =", b)
print("a + b =", overflow_result) # Wraps around to a negative number
print("Expected (120 + 10) = 130, but stored as int8:", np.int8(130))
print("\n--- Integer Overflow (int16) ---")
a16 = np.int16(32760)
b16 = np.int16(10)
print("a16 + b16 =", a16 + b16) # Correct
print("int16 max =", np.iinfo(np.int16).max)
# Show overflow explicitly
overflow16 = np.int16(np.iinfo(np.int16).max + 1)
print("Overflow: int16(max + 1) =", overflow16)
print("\n--- Integer Underflow (int8) ---")
x = np.int8(-120)
y = np.int8(-10)
underflow_result = x + y
print("x =", x)
print("y =", y)
print("x + y =", underflow_result) # Wraps around to a positive number
print("Expected (-120 - 10) = -130, but stored as int8:", np.int8(-130))
--- Integer Overflow (int8) --- a = 120 b = 10 a + b = -126
C:\Users\keith\AppData\Local\Temp\ipykernel_12860\4256458424.py:6: RuntimeWarning: overflow encountered in scalar add overflow_result = a + b
--------------------------------------------------------------------------- OverflowError Traceback (most recent call last) Cell In[15], line 10 8 print("b =", b) 9 print("a + b =", overflow_result) # Wraps around to a negative number ---> 10 print("Expected (120 + 10) = 130, but stored as int8:", np.int8(130)) 12 print("\n--- Integer Overflow (int16) ---") 13 a16 = np.int16(32760) OverflowError: Python integer 130 out of bounds for int8
import numpy as np
print("\n--- Floating-Point Overflow ---")
# This will overflow for float32 but not for float64
a = np.float32(1e38)
b = np.float32(10)
overflow_result = a * b
print("a =", a)
print("b =", b)
print("a * b =", overflow_result) # Will be inf
print("Is inf?", np.isinf(overflow_result))
print("\n--- Floating-Point Underflow ---")
# This will underflow to zero in float32
x = np.float32(1e-38)
y = np.float32(1e-10)
underflow_result = x * y
print("x =", x)
print("y =", y)
print("x * y =", underflow_result) # Will be 0.0
print("Is zero?", underflow_result == 0.0)
print("\n--- Float64 Safe Example ---")
# Compare with float64
a64 = np.float64(1e308)
b64 = np.float64(10)
print("a64 * b64 =", a64 * b64) # Will overflow
x64 = np.float64(1e-308)
y64 = np.float64(1e-10)
print("x64 * y64 =", x64 * y64) # Will underflow to subnormal or 0
--- Floating-Point Overflow --- a = 1e+38 b = 10.0 a * b = inf Is inf? True --- Floating-Point Underflow --- x = 1e-38 y = 1e-10 x * y = 0.0 Is zero? True --- Float64 Safe Example --- a64 * b64 = inf x64 * y64 = 1e-318
C:\Users\keith\AppData\Local\Temp\ipykernel_12860\1827175405.py:7: RuntimeWarning: overflow encountered in scalar multiply
overflow_result = a * b
C:\Users\keith\AppData\Local\Temp\ipykernel_12860\1827175405.py:27: RuntimeWarning: overflow encountered in scalar multiply
print("a64 * b64 =", a64 * b64) # Will overflow
import numpy as np
print("\n--- Preventing Overflow with Type Promotion ---")
# Start with int8
a = np.int8(120)
b = np.int8(10)
# Promote to int16 before addition
safe_result = np.int16(a) + np.int16(b)
print("a (int8) =", a)
print("b (int8) =", b)
print("a + b with int8 =", np.int8(a + b)) # Overflows
print("a + b with int16 =", safe_result) # Correct result 130
print("\n--- Detecting Overflow with Range Checking ---")
# Define limits
int8_max = np.iinfo(np.int8).max
int8_min = np.iinfo(np.int8).min
sum_ = int(a) + int(b)
if sum_ > int8_max or sum_ < int8_min:
print("Overflow would occur in int8: sum =", sum_)
else:
print("Safe to store in int8: sum =", sum_)
print("\n--- Use Larger Types or Unsigned Types ---")
# Promote to int32 or use uint8 if values are always positive
a_uint = np.uint8(250)
b_uint = np.uint8(10)
result_uint = a_uint + b_uint # Wraps at 255
print("uint8(250) + uint8(10) =", result_uint)
print("Expected result =", np.uint8(250 + 10)) # Will wrap to 4
# Use uint16 instead
safe_uint_result = np.uint16(a_uint) + np.uint16(b_uint)
print("Safe uint16 result =", safe_uint_result)
--- Preventing Overflow with Type Promotion --- a (int8) = 120 b (int8) = 10 a + b with int8 = -126 a + b with int16 = 130 --- Detecting Overflow with Range Checking --- Overflow would occur in int8: sum = 130 --- Use Larger Types or Unsigned Types --- uint8(250) + uint8(10) = 4
C:\Users\keith\AppData\Local\Temp\ipykernel_12860\2154227939.py:12: RuntimeWarning: overflow encountered in scalar add
print("a + b with int8 =", np.int8(a + b)) # Overflows
C:\Users\keith\AppData\Local\Temp\ipykernel_12860\2154227939.py:30: RuntimeWarning: overflow encountered in scalar add
result_uint = a_uint + b_uint # Wraps at 255
--------------------------------------------------------------------------- OverflowError Traceback (most recent call last) Cell In[16], line 32 30 result_uint = a_uint + b_uint # Wraps at 255 31 print("uint8(250) + uint8(10) =", result_uint) ---> 32 print("Expected result =", np.uint8(250 + 10)) # Will wrap to 4 34 # Use uint16 instead 35 safe_uint_result = np.uint16(a_uint) + np.uint16(b_uint) OverflowError: Python integer 260 out of bounds for uint8
Precision Errors in Dot Product¶
import numpy as np
import time
# Vector size
N = 1_000_000 # Adjust based on system memory/performance
# Generate base float64 vectors from uniform [0, 255]
A64 = (np.random.rand(N) * 255).astype(np.float64)
B64 = (np.random.rand(N) * 255).astype(np.float64)
print("\n--- float64 ---")
size_bytes = A64.nbytes + B64.nbytes
print(f"Vector size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
C64 = np.dot(A64, B64)
end_time = time.time()
print(f"Time taken: {end_time - start_time:.6f} seconds")
print("\n--- float32 ---")
A32 = A64.astype(np.float32)
B32 = B64.astype(np.float32)
size_bytes = A32.nbytes + B32.nbytes
print(f"Vector size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
C32 = np.dot(A32, B32)
end_time = time.time()
print(f"Time taken: {end_time - start_time:.6f} seconds")
print("\n--- float16 ---")
A16 = A64.astype(np.float16)
B16 = B64.astype(np.float16)
size_bytes = A16.nbytes + B16.nbytes
print(f"Vector size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
C16 = np.dot(A16, B16)
end_time = time.time()
print(f"Time taken: {end_time - start_time:.6f} seconds")
print("\n--- int32 ---")
Aint32 = A64.astype(np.int32)
Bint32 = B64.astype(np.int32)
size_bytes = Aint32.nbytes + Bint32.nbytes
print(f"Vector size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
Cint32 = np.dot(Aint32, Bint32)
end_time = time.time()
print(f"Time taken: {end_time - start_time:.6f} seconds")
print("\n--- int64 ---")
Aint64 = A64.astype(np.int64)
Bint64 = B64.astype(np.int64)
size_bytes = Aint64.nbytes + Bint64.nbytes
print(f"Vector size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
Cint64 = np.dot(Aint64, Bint64)
end_time = time.time()
print(f"Time taken: {end_time - start_time:.6f} seconds")
# Convert results to float64 for fair comparison
C32_up = np.float64(C32)
C16_up = np.float64(C16)
Cint32_up = np.float64(Cint32)
Cint64_up = np.float64(Cint64)
# Compute absolute and relative errors
print("\n--- Dot Product Absolute and Relative Errors ---")
def print_dot_error(label, C_exact, C_approx):
abs_err = np.abs(C_exact - C_approx)
rel_err = abs_err / np.abs(C_exact)
print(f"{label}: abs = {abs_err:.6e}, rel = {rel_err:.6e}")
print_dot_error("C32 vs C64", C64, C32_up)
print_dot_error("C16 vs C64", C64, C16_up)
print_dot_error("C16 vs C32", C32_up, C16_up)
print_dot_error("int32 vs C64", C64, Cint32_up)
print_dot_error("int64 vs C64", C64, Cint64_up)
--- float64 --- Vector size: 16.00 MB Time taken: 0.000976 seconds --- float32 --- Vector size: 8.00 MB Time taken: 0.000493 seconds --- float16 --- Vector size: 4.00 MB Time taken: 0.007169 seconds --- int32 --- Vector size: 8.00 MB Time taken: 0.001110 seconds --- int64 --- Vector size: 16.00 MB Time taken: 0.002229 seconds --- Dot Product Absolute and Relative Errors --- C32 vs C64: abs = 5.740566e+02, rel = 3.533729e-08 C16 vs C64: abs = inf, rel = inf C16 vs C32: abs = inf, rel = inf int32 vs C64: abs = 1.730707e+10, rel = 1.065374e+00 int64 vs C64: abs = 1.271973e+08, rel = 7.829905e-03
C:\Users\keith\AppData\Local\Temp\ipykernel_12860\3833883234.py:35: RuntimeWarning: overflow encountered in dot C16 = np.dot(A16, B16)
Precision Errors in Matrix Multiplication¶
import numpy as np
import time
# Matrix size
N = 1000 # Adjust based on system memory/performance
# Generate base float64 matrices from uniform [0, 255]
A64 = (np.random.rand(N, N) * 255).astype(np.float64)
B64 = (np.random.rand(N, N) * 255).astype(np.float64)
print("\n--- float64 ---")
size_bytes = A64.nbytes + B64.nbytes
print(f"Matrix size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
C64 = A64 @ B64
end_time = time.time()
print(f"Time taken: {end_time - start_time:.4f} seconds")
print("\n--- float32 ---")
A32 = A64.astype(np.float32)
B32 = B64.astype(np.float32)
size_bytes = A32.nbytes + B32.nbytes
print(f"Matrix size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
C32 = A32 @ B32
end_time = time.time()
print(f"Time taken: {end_time - start_time:.4f} seconds")
print("\n--- float16 ---")
A16 = A64.astype(np.float16)
B16 = B64.astype(np.float16)
size_bytes = A16.nbytes + B16.nbytes
print(f"Matrix size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
C16 = A16 @ B16
end_time = time.time()
print(f"Time taken: {end_time - start_time:.4f} seconds")
print("\n--- int32 ---")
Aint32 = A64.astype(np.int32)
Bint32 = B64.astype(np.int32)
size_bytes = Aint32.nbytes + Bint32.nbytes
print(f"Matrix size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
Cint32 = Aint32 @ Bint32
end_time = time.time()
print(f"Time taken: {end_time - start_time:.4f} seconds")
print("\n--- int64 ---")
Aint64 = A64.astype(np.int64)
Bint64 = B64.astype(np.int64)
size_bytes = Aint64.nbytes + Bint64.nbytes
print(f"Matrix size: {size_bytes / 1e6:.2f} MB")
start_time = time.time()
Cint64 = Aint64 @ Bint64
end_time = time.time()
print(f"Time taken: {end_time - start_time:.4f} seconds")
# Upsample results to float64 for fair comparison
C16_up = C16.astype(np.float64)
C32_up = C32.astype(np.float64)
Cint32_up = Cint32.astype(np.float64)
Cint64_up = Cint64.astype(np.float64)
# Compute maximum absolute differences
print("\n--- Maximum Absolute Differences ---")
print(f"max|C64 - C32| = {np.max(np.abs(C64 - C32_up)):.6e}")
print(f"max|C64 - C16| = {np.max(np.abs(C64 - C16_up)):.6e}")
print(f"max|C32 - C16| = {np.max(np.abs(C32_up - C16_up)):.6e}")
print(f"max|C64 - int32| = {np.max(np.abs(C64 - Cint32_up)):.6e}")
print(f"max|C64 - int64| = {np.max(np.abs(C64 - Cint64_up)):.6e}")
--- float64 --- Matrix size: 16.00 MB Time taken: 0.0224 seconds --- float32 --- Matrix size: 8.00 MB Time taken: 0.0077 seconds --- float16 --- Matrix size: 4.00 MB
C:\Users\keith\AppData\Local\Temp\ipykernel_12860\449388697.py:35: RuntimeWarning: overflow encountered in matmul C16 = A16 @ B16
Time taken: 7.8111 seconds --- int32 --- Matrix size: 8.00 MB Time taken: 2.1030 seconds --- int64 --- Matrix size: 16.00 MB Time taken: 2.6975 seconds --- Maximum Absolute Differences --- max|C64 - C32| = 7.997219e+00 max|C64 - C16| = inf max|C32 - C16| = inf max|C64 - int32| = 1.394282e+05 max|C64 - int64| = 1.394282e+05
Norms & Distances¶
Metric¶
Basically, boil multiple numbers (vectors, matrices, time series,...) or entire functions down to a single number
Examples
- GDP of economy
- BMI (body-mass index)
- A statistic (mean, variance, ...)
- Vector norms...
Always can be critized as discarding important info, but need to use something
Euclidean Norm¶
$$ \|x\| = \sqrt{x_1^2+x_2^2+\dotsb+x_n^2} $$"Length" of a vector (as opposed to length of the data structure, i.e. number of dimensions)
Also known as $\|x\|_2$, the "two-norm"
Example¶
$$\left\Vert \begin{pmatrix}2 \\ -1\\2 \end{pmatrix}\right\Vert = ?$$Also use a built-in function to do this
Exercise: give Euclidean norm using dot products
x = np.array([1,2,3,4])
y = np.array([5,6,7,8])
x = x.reshape(4,1) # enforce as column vectors
y = y.reshape(4,1)
print(x.T@y) # "xT*Y"
print(np.multiply(x,y))
print(x@y.T) # "xT*Y"
print(y@x.T) # "xT*Y"
print(x*y.T)
print(y*x.T)
# 1. inner product xT*y
# 2. "x*y" = np.multiply(x,y) ~ elementwise multiply
# 3.
[[70]] [[0] [0] [0] [0]] [[ 5 6 7 8] [10 12 14 16] [15 18 21 24] [20 24 28 32]] [[ 5 10 15 20] [ 6 12 18 24] [ 7 14 21 28] [ 8 16 24 32]] [[ 5 6 7 8] [10 12 14 16] [15 18 21 24] [20 24 28 32]] [[ 5 10 15 20] [ 6 12 18 24] [ 7 14 21 28] [ 8 16 24 32]]
x_sparse1 = {6:12,11:5,99:3}
for key in x_sparse1:
print(5 in x_sparse1)
False False False
Root-mean-square value (RMS)¶
$$ rms(x) = \sqrt{\frac{x_1^2+x_2^2+\dotsb+x_n^2}{n}} = \frac{1}{\sqrt{n}}\|x\| $$root(mean(square(x)))
Exercise: give Euclidean norm of sum of vectors $x+y$ in terms of dot products
(use result of previous exercise and work only with vectors)
$$||x+y||_2 = \text{dot}(?,?) + \text{dot}(?,?) + ???$$$$ \sqrt{(x_1+y_1)(x_1+y_1) + (x_2+y_2)(x_2+y_2) + ... } = $$Chebychev Inequality¶
If $x$ is a length-$n$ vector with $k$ entries satisfying $|x_i|\ge a$ for some constant $a>0$, then
$$\frac{k}{n} \le \left( \frac{rms(x)}{a}\right)^2$$Basically tells how much elements can deviate from the rms
Derive by noting that $\Vert x\Vert^2 \ge ka^2$
Motivation for other types of norms¶
Consider the list of data for different "dimensions", for one person, product, etc.
Does it make sense to compute Euclidean distance to tell us how "big" the vector is?
Norms¶
A norm is a vector "length". Often denoted generally as $\Vert \mathbf x \Vert$.
Properties
- Non-negativity: $\Vert \mathbf x \Vert \geq 0$
- Zero upon equality: $\Vert \mathbf x \Vert = 0 \iff \mathbf x = \mathbf 0$
- Absolute scalability: $\Vert \alpha \mathbf x \Vert = |\alpha| \Vert \mathbf x \Vert$ for scalar $\alpha$
- Triangle Inequality: $\Vert \mathbf x + \mathbf z\Vert \leq \Vert \mathbf x\Vert + \Vert \mathbf z\Vert$
$p$-Norms¶
For any $n$-dimensional real or complex vector. i.e. $x \in \mathbb{R}^n \text{ or } \mathbb{C}^n$
$$ \|x\|_p = \left(|x_1|^p+|x_2|^p+\dotsb+|x_n|^p\right)^{\frac{1}{p}} $$$$ \|x\|_p = \begin{pmatrix}\sum_{i=1}^n{|x_i|^p} \end{pmatrix}^{\frac{1}{p}} $$Consider the norms we have looked at. What is $p$?
Famous "norms"¶
- $\ell_2$ norm
- $\ell_1$ norm
- $\ell_\infty$ norm
- $\ell_p$ norm
- "$\ell_0$" norm
Note we often lazily write these as e.g. "L2" norm
Exercise: Wite out the $p$-norms for a vector $x$ for $p$ = 1,2,0,$\infty$... $ \|x\|_p = \begin{pmatrix}\sum_{i=1}^n{|x_i|^p} \end{pmatrix}^{\frac{1}{p}} $
Exercise: test the conditions for each
import numpy as np
v = [1,3,1,4]
for p in range(1,10):
print(p,np.power(sum(np.power(np.abs(np.array(v)),p)),1/p))
1 9.0 2 5.196152422706632 3 4.530654896083492 4 4.290915128445443 5 4.175344598847825 6 4.110988070009078 7 4.0723049678331895 8 4.048006070825583 9 4.032310478684122
Exercise¶
What are the norms of $\vec{a} = \begin{bmatrix}1\\3\\1\\-4\end{bmatrix}$ and $\vec{b} = \begin{bmatrix}2\\0\\1\\-2\end{bmatrix}$?
"Norm balls"¶
The Set $\{x | \Vert x \Vert_p \le 1 \}$, $ \|x\|_p = \begin{pmatrix}\sum_{i=1}^n{|x_i|^p} \end{pmatrix}^{\frac{1}{p}} $
Note how $\Vert x \Vert_a \le \Vert x \Vert_b$ if $a\ge b$
Norms, Inner products, and angles¶
Inner product = length squared
$$v \cdot v = v^T v = \| v \|_2^2 = \| v \|_2\| v \|_2$$Angle $\theta$ between two vectors
$$v \cdot w = v^T w = \| v \|_2\| w \|_2 \cos\theta$$Cauchy-Schwartz Inequality $v^T w \le \| v \|\| w \|$¶
Triangle Inequality¶
$$\| v + w \| \le \| v \|+\| w \|$$Derive by noting that $\| v + w \|^2 = \|v\|^2 + v^Tw + \|w\|^2$ then use Cauchy-Schwartz and complete square
$S$-norm¶
For any symmetric positive definite matrix $S$
$$\| v \|_S^2 = v^T S v$$$S$ inner product¶
$$<v,w>_S = v^T S w$$Minkowski metric¶
Proposed by Lorentz for 4-dimensional spacetime, $v = (x,y,z,t)^T$, $c$ is speed of light
$$\| v\|^2_M = x^2+y^2+z^2-ct^2$$Is it a true norm?
II. Distances¶
Euclidean distance between two vectors $a$ and $b$ in $\mathbb{R}^{n}$:¶
$$d(\mathbf a,\mathbf b) = \sqrt{\sum_{i=1}^{n}(b_i-a_i)^2} = \Vert \mathbf b - \mathbf a \Vert_2$$Again this may not make sense for various vectors. Many alternatives...
Norm versus Distance¶
What is the relationship?
Distance properties¶
A distance metic $d(\mathbf x,\mathbf y)$ must satisfy four particular conditions to be considered a metric:
- Non-negativity: $d(\mathbf x,\mathbf y) \geq 0$
- Zero upon equality: $d(\mathbf x,\mathbf y) = 0 \iff \mathbf x = \mathbf y$
- Commutativity of arguments: $d(\mathbf x,\mathbf y) = d(\mathbf y,\mathbf x)$
- Triangle Inequality: $d(\mathbf x,\mathbf z) \leq d(\mathbf x,\mathbf y) + d(\mathbf y,\mathbf z)$
Exercise¶
Write the Euclidean distance between two points entirely in terms of dot products.
What does this tell you about using dot products to compare vector similarity?
Triangle Inequality¶
$\Vert \mathbf x + \mathbf z\Vert \leq \Vert \mathbf x\Vert + \Vert \mathbf z\Vert$
$d(\mathbf a,\mathbf c) \leq d(\mathbf a,\mathbf b) + d(\mathbf b,\mathbf c)$
Application: feature distances¶
Recommender syatem based on most similar customers to movie properties
Question: what are the units of the distance here?
Application: rms prediction error¶
Time series prediction (stocks or temperature) compared to truth in retrospect
Manhattan or "Taxicab" Distance, also "Rectilinear distance"¶
Measures the relationships between points at right angles, meaning that we sum the absolute value of the difference in vector coordinates.
This metric is sensitive to rotation.
$$d_{M}(a,b) = \sum_{i=1}^{n}|b_i-a_i|$$...consider
- Non-negativity: $d(\mathbf x,\mathbf y) \geq 0$
- Zero upon equality: $d(\mathbf x,\mathbf y) = 0 \iff \mathbf x = \mathbf y$
- Commutativity of arguments: $d(\mathbf x,\mathbf y) = d(\mathbf y,\mathbf x)$
- Triangle Inequality: $d(\mathbf x,\mathbf z) \leq d(\mathbf x,\mathbf y) + d(\mathbf y,\mathbf z)$
Exercise¶
Does it fulfill the 4 conditions?
Chebyschev Distance¶
The Chebyschev distance or sometimes the $L^{\infty}$ metric, between two vectors is simply the the greatest of their differences along any coordinate dimension:
$$d_{\infty}(\mathbf a,\mathbf b) = \max_{i}{|(b_i-a_i)|}$$...consider
- Non-negativity: $d(\mathbf x,\mathbf y) \geq 0$
- Zero upon equality: $d(\mathbf x,\mathbf y) = 0 \iff \mathbf x = \mathbf y$
- Commutativity of arguments: $d(\mathbf x,\mathbf y) = d(\mathbf y,\mathbf x)$
- Triangle Inequality: $d(\mathbf x,\mathbf z) \leq d(\mathbf x,\mathbf y) + d(\mathbf y,\mathbf z)$
Cosine Distance¶
High school geometry: $\mathbf a \cdot \mathbf b = \|\mathbf a\|_2\|\mathbf b\|_2 \cos\theta$
Only depends on angle between the vectors
$$d_{\cos}(\mathbf a,\mathbf b) = 1-\frac{\mathbf a \cdot \mathbf b}{\|\mathbf a\|\|\mathbf b\|} = 1 - \cos\theta$$Not a true distance metric. which propery fails to hold? (easy to guess at based on geometry)
Exercise¶
Implement the metrics manually and compute distances between:
$ \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix}$ and $ \begin{bmatrix} 5 \\ 6 \\ 7 \\ 8 \end{bmatrix}$
import numpy as np
a = np.array([1,2,3,4])
b = np.array([5,6,7,8])
np.linalg.norm(a-b)
8.0
np.linalg.norm(a-b,np.inf)
4.0
np.sum(np.abs(a-b)) # 1 norm
16
max(np.abs(a-b)) # 1 norm
4
dcos = 1 - np.dot(a,b)/np.dot(a,a)**0.5/np.dot(b,b)**0.5
dcos
0.031136068373033732
III. Matrix Norms¶
u,s,vT = np.linalg.svd(A)
print(s)
[7.07106781e+00 2.07775413e-16]
A = np.array([[3,4],[0,0]])
print(A)
#print(np.linalg.inv(A))
print(np.linalg.cond(A))
[[3 4] [0 0]] inf
u,s,vT = np.linalg.svd(A)
print(s)
[5. 0.]
Matrix norms: motivation¶
Same goal as with vector norms: determining if matrices are small or large, or if two matrices are similar or different
Numerical mathematics example: you use 32-bit precision numbers to compute the inverse of a matrix. How different from the true inverse is your numerical calculation? Suppose you used 64-bit precision instead? (FYI the errors can be much larger than $2^{-32}$ or $2^{-64}$)
# https://www.mathworks.com/help/matlab/ref/cond.html
A1 = [[4.1, 2.8],
[9.7, 6.6]]
print(np.linalg.inv(A1))
A2 = [[4.1, 2.8 ],
[9.671, 6.608]]
print(np.linalg.inv(A2))
[[-66. 28.] [ 97. -41.]] [[ 472. -200. ] [-690.78571429 292.85714286]]
np.linalg.cond(A)
3.4032264466391812e+16
u,s,vT = np.linalg.svd(A) # A = U@S@V.T
s[0]/s[1]
3.4032264466391812e+16
np.linalg.cond([[1,0],[0,1]])
1.0
Matrix norms¶
The set of matrices can be viewed as a vector space (contains origin, closed under addition and scalar multiplication), hence we can define a norm analogously for lengths and distances
Conditions of a norm
- $\|A\| \ge 0$, only equal if $A$ is zero matrix
- $\|cA\| = |c| \|A\|$
- $\|AB\| \le \|A\|\|B\|$ ...new rule for matrix norms
Matrix norms: Froenius Norm¶
- Frobenius Norm $\Vert A \Vert_F = \sqrt{\sigma_1^2+..+\sigma_r^2} = \sqrt{\sum_{ij}A_{ij}^2}$
Exercise: Compute $\|AB\|_F$
Orthogonal matrix¶
Doesn't change L2 norm of vector or Frobenius norm matrix
$$\Vert Qx\Vert_2 = \Vert x\Vert_2$$$$\Vert QB\Vert_F = \Vert B\Vert_F$$Matrix norms from vector norms via vectorization¶
Just "vectorize" the matrix then apply a vector norm.
What do you get if you use the L2 vector norm for this?
Matrix norms from vector norms via matrix-vector product¶
Given a vector norm $\|v\|_\psi$, we can always make a matrix norm as follows
\begin{align} \|A\|_\psi &= \max_{x \ne 0}\dfrac{\Vert Ax \Vert_\psi}{\Vert x\Vert_\psi} \\ &= \max_{x \ne 0, \Vert x\Vert_\psi = 1}\Vert Ax \Vert_\psi \end{align}What does that give?¶
- $\ell_2$ $\rightarrow$ $\|A\|_2 = \sigma_1$, largest singular value of $A$
- $\ell_1$ $\rightarrow$ $\|A\|_1 = \max_{column_i}\|column_i\|_1 $, largest absolute column sum
- $\ell_\infty$ $\rightarrow$ $\|A\|_\infty = \max_{row_i}\|row_i\|_1 $, largest absolute row sum
Connections¶
$$\|A\|_\infty = \|A^T\|_1$$$$\|A\|_2^2 \le \|A\|_1 \|A\|_\infty$$Spectral and related norms¶
- Spectral Norm $\Vert A \Vert_2 = \max \dfrac{\Vert Ax \Vert}{\Vert x\Vert} = \sigma_1$, a.k.a. $\ell_2$ norm
- Frobenius Norm $\Vert A \Vert_F = \sqrt{\sigma_1^2+..+\sigma_r^2} = \sqrt{\sum_{ij}A_{ij}^2}$
- Nuclear norm $\Vert A \Vert_N = \sigma_1+..+\sigma_r$, a.k.a. trace norm
Exercise: compute norms for $I$
Spectral Radius¶
$$|\lambda|_{\max} = \max_i |\lambda_i|$$Useful metric that does not fulfill properties of norm
For every posible (true) norm, we have $\|A\| \ge |\lambda|_{\max}$
Useful fact: if $|\lambda|_{\max} <1$, $\|A^k\| \rightarrow 0$ for large $k$
Can you explain why this happens?
"Medical norm"¶
$$|||A|||_\infty = \max_{i,j}|A_{i,j}|$$Not a true matrix norm unless scaled properly also
See Strang problem I.11 #12.
Condition number¶
Quantifies how much output changes relative to a change in input.
Low condition number - good - "well conditioned"
- a small chance in inputs due to noise or numerical precision only causes a small change in output
Huge condition number - bad - "ill conditioned"
- even a tiny change in inputs can cause a drastic change in outputs. result unreliable
Condition number for solving linear system $Ax=b$¶
- View the function as $x = A^{-1}b$ where input = $b$ and output = $x$.
- Relative fractional error
Taking the max over $b$ and $\delta$ gives
\begin{align} C &= \max_{b,\delta} \dfrac{\|A^{-1}\delta\|}{\|A^{-1}b\|}\dfrac{\|b\|}{\|\delta\|} \\ &= \left(\max_{\delta}\dfrac{\|A^{-1}\delta\|}{\|\delta\|}\right) \left(\max_{b}\dfrac{\|b\|}{\|A^{-1}b\|} \right) \\ &= \|A^{-1}\|\|A\| \end{align}First term is just definition of matrix norm. For second term, plug in $b = Ax$ and get norm definition also, noting that max will be the same optimizing over $x$ or over $Ax$ since $A$ is full-rank.
"Matrix Condition number"¶
Usually means condition number with L2 norm
Recall $\|A\|_2 = \sigma_\max(A)$
$$C \equiv \|A^{-1}\|_2\|A\|_2 = \dfrac{\sigma_\max(A)}{\sigma_\min(A)}$$Suppose a matrix is singular? (what does this mean?)
- what are its singular values?
- what is its condition number?
# https://www.mathworks.com/help/matlab/ref/cond.html
A1 = [[4.1, 2.8],
[9.7, 6.6]]
print(np.linalg.inv(A1))
A2 = [[4.1, 2.8 ],
[9.671, 6.608]]
print(np.linalg.inv(A2))
[[-66. 28.] [ 97. -41.]] [[ 472. -200. ] [-690.78571429 292.85714286]]
print(np.linalg.cond(A1),np.linalg.cond(A2))
1622.9993838565106 11560.27884206701
v1 = np.array([1.,2.])
A = np.array([v1,2.*v1])
print(A)
print('rank =',np.linalg.matrix_rank(A))
[[1. 2.] [2. 4.]] rank = 1
An = A + 1e-10*np.random.rand(2,2)
print('rank =',np.linalg.matrix_rank(An))
print('cond =',np.linalg.cond(An))
rank = 2 cond = 282960980020.97266
Exercise¶
Test the prior matrices when solving a linear system $Ax=b$ and adding a small perturbation to $b$
Change the matrix to make it weakly conditioned and repeat,
# Compute maximum absolute differences
print("\n--- Maximum Absolute Differences ---")
print(f"max|C64 - C32| = {np.max(np.abs(C64 - C32_up)):.6e}")
print(f"max|C64 - C16| = {np.max(np.abs(C64 - C16_up)):.6e}")
print(f"max|C32 - C16| = {np.max(np.abs(C32_up - C16_up)):.6e}")
print(f"max|C64 - int32| = {np.max(np.abs(C64 - Cint32_up)):.6e}")
print(f"max|C64 - int64| = {np.max(np.abs(C64 - Cint64_up)):.6e}")
# Compute Frobenius norm absolute and relative errors
print("\n--- Frobenius Norm Errors ---")
def print_frobenius_error(label, C_exact, C_approx):
abs_err = np.linalg.norm(C_exact - C_approx, ord='fro')
rel_err = abs_err / np.linalg.norm(C_exact, ord='fro')
print(f"{label}: abs = {abs_err:.6e}, rel = {rel_err:.6e}")
print_frobenius_error("C32 vs C64", C64, C32_up)
print_frobenius_error("C16 vs C64", C64, C16_up)
print_frobenius_error("C16 vs C32", C32_up, C16_up)
print_frobenius_error("int32 vs C64", C64, Cint32_up)
print_frobenius_error("int64 vs C64", C64, Cint64_up)
--- Maximum Absolute Differences --- max|C64 - C32| = 7.997219e+00 max|C64 - C16| = inf max|C32 - C16| = inf max|C64 - int32| = 1.394282e+05 max|C64 - int64| = 1.394282e+05 --- Frobenius Norm Errors --- C32 vs C64: abs = 1.599064e+03, rel = 9.823474e-08 C16 vs C64: abs = inf, rel = inf C16 vs C32: abs = inf, rel = inf int32 vs C64: abs = 1.273148e+08, rel = 7.821287e-03 int64 vs C64: abs = 1.273148e+08, rel = 7.821287e-03
IV. Statistics with Python¶
Random sampling¶
- rand() - uniform random numbers in [0,1]
- randn() - standard normal random numbers (zero mean, unit variance)
from numpy import random
printcols(dir(random))
BitGenerator __package__ default_rng noncentral_chisquare set_state Generator __path__ dirichlet noncentral_f shuffle MT19937 __spec__ exponential normal standard_cauchy PCG64 _bounded_integers f pareto standard_exponential PCG64DXSM _common gamma permutation standard_gamma Philox _generator geometric poisson standard_normal RandomState _mt19937 get_state power standard_t SFC64 _pcg64 gumbel rand test SeedSequence _philox hypergeometric randint triangular __RandomState_ctor _pickle laplace randn uniform __all__ _sfc64 logistic random vonmises __builtins__ beta lognormal random_integers wald __cached__ binomial logseries random_sample weibull __doc__ bit_generator mtrand ranf zipf __file__ bytes multinomial rayleigh __loader__ chisquare multivariate_normal sample __name__ choice negative_binomial seed
The Normal Distribution¶
- Also known as Gaussian distribution
import numpy as np
from matplotlib import pyplot as plt
def univariate_normal(x, mean, var):
return ((1. / np.sqrt(2 * np.pi * var)) * np.exp(-(x - mean)**2 / (2 * var)))
x = np.linspace(-5,5,1000)
plt.plot(x,univariate_normal(x,1,2));
plt.show()
Exercise¶
Generate normal random samples with mean of 1 and variance of 1
Plot histogram versus the theoretical distribution
Try different numbers of samples
Statistics¶
Consider the relation between norms and simple statistical quantities
\begin{align} \text{Population mean} &= \mu = \frac{\sum_{i=1}^N x_i}{N} \\ \text{Sample mean} &= \bar{x} = \frac{\sum_{i=1}^n x_i}{n} \\ \text{Population variance} &= \sigma^2 = \frac{\sum_{i=1}^N (x_i - \mu)^2}{N} \\ \text{Sample variance} &= s^2 = \frac{\sum_{i=1}^n (x_i - \bar{x})^2}{n - 1} = \frac{\sum_{i=1}^n x_i^2 - \frac{1}{n}(\sum_{i=1}^n x_i)^2}{n - 1} \\ \text{Standard deviation} &= \sqrt{\text{Variance}} \end{align}Exercise¶
Load a dataset from sklearn and compute mean and variance of a column
Standardize the column
Now compute mean and variance of the result
The Standard Normal Distribution $Z$¶
- If $X$ is a normal r.v. with $E(X) = \mu$ and $V(X) = \sigma^2$
- $Z$ is a normal r.v. with $E(X) = 0$ and $V(X) = 1$
Standardizing data has two tasks¶
$$ Z = \frac{X-\mu}{\sigma} $$- Remove the mean
- scale by the standard deviation
Lab: Standardizing data¶
Standardize the columns of the Iris dataset using linear algebra.
Test it worked by computing the mean and norm of each column.
Multivariate Gaussian (for $n$ dimensions)¶
$$ f(\mathbf x) = \frac{1}{ \sqrt{2 \pi^n |\boldsymbol\Sigma|}} \exp \left(- \frac{1}{2} (\mathbf x - \boldsymbol \mu)^T \boldsymbol\Sigma^{-1} (\mathbf x - \boldsymbol \mu) \right) \text{, for } \mathbf x \in R^n $$- Mean vector as centroid of distribution
- Covariance matrix describes spread - correlations between variables $\Sigma_{ij} = S_{\mathbf x_i \mathbf x_j}$
def multivariate_normal(x, n, mean, cov):
return (1./(np.sqrt((2*np.pi)**n * np.linalg.det(cov))) * np.exp(-1/2*(x - mean).T@np.linalg.inv(cov)@(x - mean)))
mean = np.array([35,70])
cov = 100*np.array([[1,.5],[.5,1]])
pic = np.zeros((100,100))
for x1 in np.arange(0,100):
for x2 in np.arange(0,100):
x = [x1,x2]
pic[x1,x2] = multivariate_normal(x, 2, mean, cov)
plt.contour(pic);
Tricky Exercise¶
Consider a $m \times 3$ matrix $\mathbf A$ with columns $\mathbf x, \mathbf y$, and $\mathbf z$ (each a vector of data).
How would you efficiently standardize the three columns to make $\bar{\mathbf A}$.
What are the elements of $\bar{\mathbf A}^T \bar{\mathbf A}$?
Exercise:¶
Assume you have a vector containing samples. Write the following in terms of norms and dot products:
- mean
- variance
- Correlation coefficient
- Covariance
So what does this tell you about comparing things using distances versus dot products versus statistics?
V. Clustering¶
Dimensionality Reduction¶
Describing data with less information
How are PCA and Clustering leading to less information? Consider a dataset.
What are the benefits?
Clustering References¶
https://en.wikipedia.org/wiki/K-means_clustering
https://en.wikipedia.org/wiki/K-means%2B%2B - kmeans++
https://www.youtube.com/watch?v=IuRb3y8qKX4 - video with visualization of training progress
https://www.youtube.com/watch?v=cWSnFaSjgBU - more on visualization
http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html
Marketing Motivation¶
You want to make a certain number of products for a large population of customers. We know a number of features describing each of the customers, and are able to target a product to a certain customer profile.
E.g.: different customers prefer cars which are fast, or get good gas mileage, or are cheap, or are luxurious and presigious (i.e., expensive), or are big, or are little, etc. And in various combinations.
The customers vary widely so we wish to make k products which cover as closely as possible what as many as possible customers would like.
The better you can do this, the better your products will sell.
K-means clustering¶
- K clusters, centers are the means of cluster members.
- Each sample belongs to the cluster who's mean is nearest.
- Iteratively recompute membership and means (greedy optimization)
K-means algorithm¶
- Assign a number from 1 to K to each of N data points randomly
- While cluster assignments keep changing:
- For each of K clusters:
- Calculate cluster centroid
- For each of N points:
- Assign point to centroid it is closest to
- For each of K clusters:
Consider the possible ways to vary details in this method? How many clusters? How calculate distance?
k-Means as Greedy Optimization¶
The sum of all distances between members of a cluster gives a total length (a squared length can be viewed as a measure of "energy").
Within-Cluster-Variation: $WCV(C_k) = \dfrac{1}{|C_k|}\sum_{i, j \in C_k}d(x_{i}-x_{j})$, where $d(x_{i}-x_{j})$ is a distance metric of your choice.
k-means tries to minimize the net lengths over all clusters.
If you plotted this for each iteration, what would the plot look like?
Lab¶
Implement k-means clustering algorithm using scikit on IRIS, MNIST datasets.
Compute the net WCV for various options.
Look at result and WCV with different initializations.
Try varying the metric used in the clustering.
Vary k and plot WCV versus k.