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 0x7f43a5708450>
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 0x7f43a570a650>

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.

1. We should only allow Rational to be constructed from integers

2. 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)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[8], line 1
----> 1 Rational(1.0, 2.0)

Cell In[7], line 16, in Rational.__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)
Cell In[9], line 1
----> 1 Rational(1,0)

Cell In[7], line 14, in Rational.__init__(self, p, q)
      9 """
     10 initialize a Rational object with numerator p, denominator q
     11 """
     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

1. __float__ to return the floating point representation of a Rational

2. __mul__ to overload the multiplication operator x * y

3. __truediv__ to overload the division operator x / y

4. __sub__ to overload the subtraction operator x - y

5. __neg__ to overload the negation operator -x

6. __abs__ to return the absolute value abs(x)

Implement some checks to make everything works.

Now, you can do rational arithmetic!

# Your code here

Show code cell content Hide code cell content

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)
Cell In[13], line 3
      1 A = LinearMap()
      2 x = np.random.rand(2)
----> 3 A * x

Cell In[12], line 46, in LinearMap.__mul__(self, other)
     43     return ProdMap(other, self)
     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')

Cell In[12], line 68, in LinearMap.matmul(self, x)
     66     return self.data @ x # numpy matrix-vector multiplication
     67 else:
---> 68     raise NotImplementedError("matmul")

NotImplementedError: matmul

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

LinearMap(<class 'numpy.ndarray'>)

array([0.60738397, 0.14364453])

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()

(-35.03831609111494, -35.038316091114865)

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.01990770300000122 sec.
time elapsed: 0.0008204000000091582 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

1. Implement a derived class of LinearMap for the identity map \(I: x\to x\), called IdentityMap

2. 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 = 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

Show code cell content Hide code cell content

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)

Show code cell content Hide code cell content

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]
    

Show code cell content Hide code cell content

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.0

Show code cell content Hide code cell content

x0 = np.random.randn(n)
RayleighQuotient(A, x0)

1.0042892412341002

Show code cell content Hide code cell content

x1, lam = PowerMethod(A, x0)