Object-Oriented Programming - Scientific Computing with Python

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}")

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)

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(f'x = {x}')
x

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)

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)

We should also check that the inputs to our Rational class are valid.

We should only allow Rational to be constructed from integers
The 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)

Rational(1,0)

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 operator x * y
__truediv__ to overload the division operator x / y
__sub__ to overload the subtraction operator x - y
__neg__ to overload the negation operator -x
__abs__ to return the absolute value abs(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):

If $f: \mathbb R^m \to \mathbb R^n$ , $g : \mathbb R^m \to \mathbb R^n$ , then $f + g : \mathbb R^m \to \mathbb R^n$ is also a linear map.
If $f : \mathbb R^m \to \mathbb R^n$ , $h : \mathbb R^n \to \mathbb R^p$ , then $f \circ h : \mathbb R^m \to \mathbb R^p$ is also a linear map.

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, that takes in a vector and returns the output. This will only work after we pass in a data parameter.

A = LinearMap()
x = np.random.rand(2)
A * x

A = LinearMap(np.random.rand(2,2))
print(A)

A * x

The @property decorator makes a method accessible as a property

A.shape

To finish the LinearMap functionality, we need to define some derived classes (child classes) ProdMap and SumMap. Since the sum and product of linear maps are still linear maps, we want to inherit from our LinearMap parent class.

class ProdMap(LinearMap):
    """
    Product of Linear Maps (A * B)
    """
    def __init__(self, A, B): # Notice it requires two inputs for initialization
        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): # The matmul function is overriden
        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)

M = LinearMap(np.ones((2,2))) * 2
N = LinearMap(np.eye(2))

A = M + N
x = np.array([1,0])

A * x

Note that we can do arithmetic on SumMap and ProdMap

B = A * A
B * x

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
M = np.random.randn(n, k)
N = np.random.randn(k, n)

dense_matrix = M @ N # dense matrix product
print(type(dense_matrix))

linear_map = LinearMap(M) * LinearMap(N)
print(linear_map)

dense_matrix.trace(), linear_map.trace()

import time

x = np.random.randn(n)
t0 = time.monotonic()
for i in range(10):
    y = dense_matrix @ x
t1 = time.monotonic()
print("time elapsed: {} sec.".format(t1 - t0))

t0 = time.monotonic()
for i in range(10):
    y = linear_map * x
t1 = time.monotonic()
print("time elapsed: {} sec.".format(t1 - t0))

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 LinearMap for the identity map $I: x\to x$ , called IdentityMap
Implement a derived class of LinearMap called SymmetricMap which 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 = \lambda * x$ . 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}

(1)

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}

(2)

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]

A.trace()