RLS.py

""" Code that executes a recursive least squares algorithm """

import numpy as np
import os

class RecursiveLeastSquaresError(Exception):
    pass

class RecursiveLeastSquares:
    def __init__(self, num_iterations, num_dimensions, plotting=False):
        self.num_iterations = num_iterations
        if self.num_iterations < 0 or self.num_iterations > 10:
            raise RecursiveLeastSquaresError(f"Number of iterations {self.num_iterations} if out of the supported range. Please enter a number 1-9.")
        self.num_dimensions = num_dimensions
        if self.num_dimensions > 60:
            raise RecursiveLeastSquaresError("Number of dimensions is too high for this example. Choose something less than 60.")
        self.plotting = plotting

        self.x = np.random.rand(self.num_dimensions, 1)
        self.y = np.random.rand(self.num_dimensions, 1)
        self.A = np.column_stack((np.ones(len(self.x)), self.x))
        self.b = self.y.reshape(-1, 1)

    def least_squares(self, confirm_errors=False):
        """ Implements the normal equations from my previous project (https://github.com/zaynpatel/TAAS/blob/main/Deliverable_1/Least_Squares.ipynb) """
        A_T_A = np.linalg.inv(self.A.T @ self.A)
        A_T_b = self.A.T @ self.b
        x_hat = A_T_A @ A_T_b

        if confirm_errors:
            projection = self.A @ x_hat
            error = self.b - projection
            residuals = sum([res[0] * res[0] for res in error])
            print(f"The residual error is: {residuals}")
            assert self.b.all() == projection.all() + error.all()
            np_soln = np.linalg.lstsq(self.A, self.b, rcond=None)
            assert x_hat.all() == np_soln[0].all()
            return x_hat

        return x_hat

    # function that runs recursive least squares and nicely prints the new slope and intercept values
    def recursive_least_squares(self):
        np.random.seed(4)
        if self.num_iterations == 0:
            x_hat = self.least_squares()
            return x_hat

        for _ in range(self.num_iterations):
            x_new_row = np.random.rand(1, 2)
            b_new = np.random.rand(1, 1)
            x_new_equation = self.sherman_morrison_formula(x_new_row) @ ((self.A.T @ self.b) + (x_new_row.T * b_new))
            print(f"New equation for rank 1 change to A transpose A: {x_new_equation}")
            self.A = np.append(self.A, x_new_row, axis=0)
            self.b = np.append(self.b, b_new, axis=0)
            np_rls_answer = np.linalg.lstsq(self.A, self.b, rcond=None)
            print(f"Final RLS answer: {np_rls_answer[0]}")
            assert np_rls_answer[0].all() == x_new_equation.all()

        return

    def sherman_morrison_formula(self, x_new_row):
        """ Implements a version of the Sherman-Morrison formula (https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula) """
        ATA_inv = np.linalg.inv(self.A.T @ self.A)
        y = ATA_inv @ np.transpose(x_new_row)
        c = 1/(1 + (x_new_row) @ y)
        return ATA_inv - (c * ATA_inv @ np.transpose(x_new_row) @ x_new_row @ ATA_inv)


    # function for plotting (I want capability to plot multiple figures in one place)
    def plot(self):
        pass