Object-Oriented Programming
Contents
Object-Oriented Programming¶
Object-oriented programming is a style of programming that bundles data into objects and defines things that might be done to objects using methods.
Programs written in an object-oriented style may look very different from programs written in a procedural or funtctional style. If you’re coming to Python from Matlab, Fortran, or C, this can take some time to get used to.
Classes¶
A Python class is a recipe for creating objects of a certain type. An object created from a class is an instantiation of the class.
Typically, there are operations you may like to perform on objects in a certain class. You can define these operations as functions declared within the class, called methods.
Classes typically hold some sort of data. The __init__() method defines how this data is stored. To allow an object to refer to its own data, we typically use self
class Rational(object):
"""
A class for rational numbers p / q
"""
def __init__(self, p, q):
"""
initialize a Rational object with numerator p, denominator q
"""
self.p = p # store p as a parameter
self.q = q # store q as a parameter
x = Rational(1, 2)
print(x)
print(f"x = {x.p} / {x.q}")
<__main__.Rational object at 0x7fdb8cebb970>
x = 1 / 2
Clearly, we might like to have a method that turns a Rational into a nicely formatted string. We might accomplish by adding a method tostring
class Rational(object):
"""
A class for rational numbers p / q
"""
def __init__(self, p, q):
"""
initialize a Rational object with numerator p, denominator q
"""
self.p = p # store p as a parameter
self.q = q # store q as a parameter
def tostring(self):
return f"{self.p} / {self.q}"
x = Rational(1, 2)
print(x.tostring())
print(x)
1 / 2
<__main__.Rational object at 0x7fdb8cea3790>
However, the tostring method doesn’t play nicely with things like print. The solution is to use one of Python’s dunder (“magic”) methods. These are special methods, always surrounded by a double underscore (“dunder” = “Double UNDERscore”) like __init__, which Python knows to look for in certain situations.
The magic method to turn an object into a string is the __str__() method.
The __repr__() method is the official representation of the object.
class Rational(object):
"""
A class for rational numbers p / q
"""
def __init__(self, p, q):
"""
initialize a Rational object with numerator p, denominator q
"""
self.p = p # store p as a parameter
self.q = q # store q as a parameter
def __str__(self):
return f"{self.p} / {self.q}"
def __repr__(self):
return f"Rational({self.p}, {self.q})"
x = Rational(1, 2)
print('x = {}'.format(x))
x
x = 1 / 2
Rational(1, 2)
You can find documentation on the different magic methods you might use in the Python documentation here.
One of the things we would often like to do in scientific computing is define types of mathematical objects and do operations on them in a natural way.
class Rational(object):
"""
A class for rational numbers p / q
"""
def __init__(self, p, q):
"""
initialize a Rational object with numerator p, denominator q
"""
self.p = p # store p as a parameter
self.q = q # store q as a parameter
def __str__(self):
return f"{self.p} / {self.q}"
def __repr__(self):
return f"Rational({self.p}, {self.q})"
def __add__(self, other):
"""
Rational addition, self + other
"""
return Rational(self.p * other.q + self.q * other.p, self.q * other.q)
x = Rational(1,2)
y = Rational(1,4)
z = x + y
print(z)
6 / 8
we might want to add a new method that simplifies the fraction by removing any common divisors from the numerator and denominator, and ensuring that the denominator is positive.
from math import gcd
class Rational(object):
"""
A class for rational numbers p / q
"""
def __init__(self, p, q):
"""
initialize a Rational object with numerator p, denominator q
"""
self.p = p # store p as a parameter
self.q = q # store q as a parameter
self.simplify()
def __str__(self):
return f"{self.p} / {self.q}"
def __repr__(self):
return f"Rational({self.p}, {self.q})"
def __add__(self, other):
"""
Rational addition
"""
return Rational(self.p * other.q + self.q * other.p, self.q * other.q)
def simplify(self):
"""
simplify fraction to lowest terms, and ensure denominator is positive
"""
g = gcd(self.p, self.q)
sgn = 1 if self.q > 0 else -1 # ternary expression in Python
self.p = sgn * self.p // g
self.q = sgn * self.q // g
x = Rational(1,2)
y = Rational(1,4)
z = x + y
print(z)
w = Rational(1,-2)
print(w)
3 / 4
-1 / 2
We should also check that the inputs to our Rational class are valid.
We should only allow
Rationalto be constructed from integersThe denominator should not be zero
we can use the isinstance function to determine if a variable is of a certain class.
from math import gcd
class Rational(object):
"""
A class for rational numbers p / q
"""
def __init__(self, p, q):
"""
initialize a Rational object with numerator p, denominator q
"""
if q == 0:
raise ValueError('Denominator must not be zero')
if not isinstance(p, int):
raise ValueError('Numerator must be an integer')
if not isinstance(q, int):
raise ValueError('Denominator must be an integer')
self.p = p # store p as a parameter
self.q = q # store q as a parameter
self.simplify()
def __str__(self):
return f"{self.p} / {self.q}"
def __repr__(self):
return f"Rational({self.p}, {self.q})"
def __add__(self, other):
"""
Rational addition
"""
return Rational(self.p * other.q + self.q * other.p, self.q * other.q)
def simplify(self):
"""
simplify fraction to lowest terms, and ensure denominator is positive
"""
g = gcd(self.p, self.q)
sgn = 1 if self.q > 0 else -1 # ternary expression in Python
self.p = sgn * self.p // g
self.q = sgn * self.q // g
Rational(1.0, 2.0)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-8-41c18925afd4> in <module>
----> 1 Rational(1.0, 2.0)
<ipython-input-7-6cac91fe9a82> in __init__(self, p, q)
14 raise ValueError('Denominator must not be zero')
15 if not isinstance(p, int):
---> 16 raise ValueError('Numerator must be an integer')
17 if not isinstance(q, int):
18 raise ValueError('Denominator must be an integer')
ValueError: Numerator must be an integer
Rational(1,0)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-9-e6a8013a20b2> in <module>
----> 1 Rational(1,0)
<ipython-input-7-6cac91fe9a82> in __init__(self, p, q)
12
13 if q == 0:
---> 14 raise ValueError('Denominator must not be zero')
15 if not isinstance(p, int):
16 raise ValueError('Numerator must be an integer')
ValueError: Denominator must not be zero
Exercise¶
Implement the following magic methods for the Rational class
__float__to return the floating point representation of a Rational__mul__to overload the multiplication operatorx * y__truediv__to overload the division operatorx / y__sub__to overload the subtraction operatorx - y__neg__to overload the negation operator-x__abs__to return the absolute valueabs(x)
Implement some checks to make everything works.
Now, you can do rational arithmetic!
# Your code here
class Rational(object):
"""
A class for rational numbers p / q
"""
def __init__(self, p, q):
"""
initialize a Rational object with numerator p, denominator q
"""
if q == 0:
raise ValueError('Denominator must not be zero')
if not isinstance(p, int):
raise ValueError('Numerator must be an integer')
if not isinstance(q, int):
raise ValueError('Denominator must be an integer')
self.p = p # store p as a parameter
self.q = q # store q as a parameter
self._simplify()
def __str__(self):
return f"{self.p} / {self.q}"
def __repr__(self):
return f"Rational({self.p}, {self.q})"
def __add__(self, other):
"""
Rational addition
"""
return Rational(self.p * other.q + self.q * other.p, self.q * other.q)
def __float__(self):
return self.p / self.q
def __mul__(self, other):
return Rational(self.p * other.p, self.q * other.q)
def __truediv__(self, other):
return Rational(self.p * other.q, self.q * other.p)
def __sub__(self, other):
return Rational(self.p * other.q - self.q * other.p, self.q * other.q)
def __neg__(self):
return Rational(-self.p, self.q)
def __abs__(self):
return Rational(abs(self.p), self.q)
def _simplify(self):
"""
simplify fraction to lowest terms, and ensure denominator is positive
"""
g = gcd(self.p, self.q)
sgn = 1 if self.q > 0 else -1 # ternary expression in Python
self.p = sgn * self.p // g
self.q = sgn * self.q // g
Class Inheritance¶
Sometimes, multiple classes are related even though they have different data. In mathematics, we might say that the objects are in the same group/ring etc.
One example is a class of linear maps \(f: \mathbb{R}^m \to \mathbb{R}^n\). Linear maps can be added, scaled, or composed (multiplied).
Linear maps can generally be defined using a dense matrix, but that isn’t always the most efficient representation.
We’ll define a base class (also called a parent class) which will define how linear maps behave.
import numpy as np
import numbers
class LinearMap(object):
"""
A Linear Map from R^m -> R^n
"""
def __init__(self, data=None):
self.data = data
@property
def shape(self):
"""
return the shape of the Linear Map
"""
if self.data is not None and type(self.data) is np.ndarray:
return self.data.shape # assume data has a shape
else:
return None, None # indeterminate shape
def __repr__(self):
return f"LinearMap({type(self.data)})"
def __add__(self, other):
"""
Addition of linear maps
"""
if isinstance(other, LinearMap) or isinstance(other, numbers.Number):
return SumMap(self, other)
else:
raise ValueError('other type must be LinearMap!')
def __mul__(self, other):
"""
Multiplication of linear maps
"""
if isinstance(other, LinearMap):
return ProdMap(self, other)
elif isinstance(other, numbers.Number):
# other is a scalar
return ProdMap(other, self)
elif isinstance(other, np.ndarray):
# other is a vector
return self.matmul(other)
else:
raise ValueError('unsupported type for multiplication')
def __rmul__(self, other):
"""
This is called if the other type is multiplied on the left
other * self
"""
if isinstance(other, numbers.Number):
# other is a scalar
return ProdMap(other, self)
else:
raise ValueError('unsupported type for multiplication')
def matmul(self, x):
"""
Implements action of the LinearMap on a vector
"""
if self.data is not None:
return self.data @ x # numpy matrix-vector multiplication
else:
raise NotImplementedError("matmul")
def trace(self):
"""
Compute the trace of LinearMap, assuming shape is determinate
"""
if self.shape[0] is not None and self.shape[1] is not None:
if self.shape[0] == self.shape[1]:
n = self.shape[0]
s = 0 # running sum
for i in range(n):
ei = np.zeros(n)
ei[i] = 1
s = s + np.dot(ei, self.matmul(ei))
return s
else:
raise RuntimeError('LinearMap is not square!')
else:
raise RuntimeError('LinearMap has indeterminate shape')
x = LinearMap()
You’ll see above that we use several classes above: ProdMap, and SumMap. We also define a matmul function, but we haven’t implemented it, unless we pass in a data parameter.
A = LinearMap()
x = np.random.rand(2)
A * x
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
<ipython-input-13-366473083d45> in <module>
1 A = LinearMap()
2 x = np.random.rand(2)
----> 3 A * x
<ipython-input-12-db0b9c4e454d> in __mul__(self, other)
44 elif isinstance(other, np.ndarray):
45 # other is a vector
---> 46 return self.matmul(other)
47 else:
48 raise ValueError('unsupported type for multiplication')
<ipython-input-12-db0b9c4e454d> in matmul(self, x)
66 return self.data @ x # numpy matrix-vector multiplication
67 else:
---> 68 raise NotImplementedError("matmul")
69
70 def trace(self):
NotImplementedError: matmul
A = LinearMap(np.random.rand(2,2))
print(A)
A * x
LinearMap(<class 'numpy.ndarray'>)
array([0.19703971, 0.07614441])
The @property decorator makes a method accessible as a property
A.shape
(2, 2)
What we need to do is define some derived classes (child classes) that inherit from our LinearMap parent class.
class ProdMap(LinearMap):
"""
Product of Linear Maps (A * B)
"""
def __init__(self, A, B):
self.A = A
self.B = B
@property
def shape(self):
return self.A.shape[0], self.B.shape[1]
def __repr__(self):
return f"ProdMap({self.A}, {self.B})"
def matmul(self, x):
return self.A * (self.B * x)
class SumMap(LinearMap):
"""
Sum of Linear Maps (A + B)
"""
def __init__(self, A, B):
self.A = A
self.B = B
@property
def shape(self):
return self.A.shape[0] or self.B.shape[0], self.A.shape[1] or self.B.shape[1]
def __repr__(self):
return f"SumMap({self.A}, {self.B})"
def matmul(self, x):
return (self.A * x) + (self.B * x)
A = SumMap(3,2)
x = np.ones((2,1))
A * x
array([[5.],
[5.]])
Note that we can do arithmetic on SumMap and ProdMap
B = A * A
print(B)
B * x
ProdMap(SumMap(3, 2), SumMap(3, 2))
array([[25.],
[25.]])
Because SumMap and ProdMap are derived classes of LinearMap, they inherit all the methods of LinearMap by default. This means that we can automatically multiply two SumMap objects without needing to define __mul__ in the derived class.
Note, that if we define a method in both the base class and derived class that the derived class method is called. This happens for the __init__, __repr__ and matmul methods above.
One of the things we get with our LinearMap class is a cheap way to construct low-rank matrices
n = 2000
k = 10
A = np.random.randn(n, k)
B = np.random.randn(k, n)
AB = A @ B # dense product
AB2 = LinearMap(A) * LinearMap(B)
AB2
ProdMap(LinearMap(<class 'numpy.ndarray'>), LinearMap(<class 'numpy.ndarray'>))
AB2.trace(), AB.trace()
(-124.50961244079777, -124.50961244079781)
import time
x = np.random.randn(n)
t0 = time.monotonic()
for i in range(10):
y = AB @ x
t1 = time.monotonic()
print("time elapsed: {} sec.".format(t1 - t0))
t0 = time.monotonic()
for i in range(10):
y = AB2 * x
t1 = time.monotonic()
print("time elapsed: {} sec.".format(t1 - t0))
time elapsed: 0.02155948698054999 sec.
time elapsed: 0.0008220209856517613 sec.
Why is the ProdMap faster? One way to reason about this is that it only stores about 1/100th of the data of the dense matrix: \((2000 \times 10 \times 2) / (2000 \times 2000) = 1/100\). Because matrix-vector multiplication requires us to look at every element of the matrix, it should take about 1/100th of the time to loop over the data in ProdMap (there’s more to consider, but this is an ok rule of thumb).
Exercise¶
Part 1 - More Derived Classes¶
Implement a derived class of
LinearMapfor the identity map \(I: x\to x\), calledIdentityMapImplement a derived class of
LinearMapcalledSymmetricMapwhich has as data a matrix \(A\), and implements the linear map \(A * A^T\)
Part 2 - Power method¶
Recall an eigenvector of a linear map \(A\) is a vector \(x\) so that \(A x = x * \lambda\). The scalar \(\lambda\) is the eigenvalue, and the pair \((\lambda, x)\) is called the eigenpair.
Now we’ll implement power method to find the largest eigenpair of a symmetric (hermetian) matrix \(A\). In pseudo-code, the algorithm is
Inputs:
A: a symmetric n x n matrix
x: a vector of length n
while not converged:
x = A * x
x = x / ||x||_2
The vector x will converge to the eigenvector with largest magnitude eigenvalue, and the eigenvalue can be computed as the Rayleigh quotient
$\(R(A, x) = \frac{x^T A x}{x^T x}\)$
If \(x_k\) is the value of \(x\) at iteration \(k\), we’ll say the algorithm has converged if
$\(|R(A, x_k) - R(A, x_{k+1}| < \mathsf{tol}\)$
where tol is some tolerance.
Write a function that will use the power method to compute the largest eigenpair. You should be able to call the function as
x, lam = PowerMethod(A, x0, tol=1e-8) # default parameters provided
where x0 is the inital vector used in the iteration.
You may want to consider writing a helper function to compute the Rayleigh quotient.
Test your function on a matrix \(A = I + x x^T\) where \(x\) is a randomly generated unit vector. How does the top eigenvector relate to \(x\)? What is the top eigenvalue? You can answer this either with your knowledge of linear algebra, or experimentally.
# Your code here
class IdentityMap(LinearMap):
"""
The identity map I : x -> x
"""
def __init__(self):
self.data = None
def __repr__(self):
return "IdentityMap()"
def matmul(self, x):
return x
class SymmetricMap(LinearMap):
"""
A map A * A.T()
"""
def __init__(self, A):
self.A = A
@property
def shape(self):
return self.A.shape[0], self.A.shape[0]
def __repr__(self):
return "SymmetricMap({})".format(type(A))
def matmul(self, x):
return self.A @ (self.A.T @ x)
def RayleighQuotient(A, x):
"""
computes Rayleigh quotient
R(A, x) = x^T A x / x^T x
"""
return x.dot(A * x) / x.dot(x)
def PowerMethod(A, x0, tol=1e-8):
"""
compute top eigenpair of A using power method. Start with vector x0
"""
x = x0
r = [RayleighQuotient(A, x)]
while True:
x = A * x
x = x / np.linalg.norm(x)
r.append(RayleighQuotient(A, x))
if abs(r[-1] - r[-2]) < tol:
break
return x, r[-1]
n = 100
x = np.random.randn(n)
x = x.reshape(n, -1)
x = x / np.linalg.norm(x)
A = IdentityMap() + SymmetricMap(x)
A.trace()
101.00000000000001
x0 = np.random.randn(n)
RayleighQuotient(A, x0)
1.0072034638701906
x1, lam = PowerMethod(A, x0)