HM1_Aufgabenserie8/Schenk_Brandenberger_S8_Aufg2.py

# -*- coding: utf-8 -*-
"""
Created on Sat Nov  7 13:26:09 2020

Höhere Mathematik 1, Serie 8, Gerüst für Aufgabe 2

Description: calculates the QR factorization of A so that A = QR
Input Parameters: A: array, n*n matrix
Output Parameters: Q : n*n orthogonal matrix
                   R : n*n upper right triangular matrix            
Remarks: none
Example: A = np.array([[1,2,-1],[4,-2,6],[3,1,0]]) 
        [Q,R]=Serie8_Aufg2(A)

@author: knaa
"""
import numpy as np
import timeit

# Aufgabe 2a
def Serie8_Aufg2(A):
    
    A = np.copy(A)                       #necessary to prevent changes in the original matrix A_in
    A = A.astype('float64')              #change to float
    
    n = np.shape(A)[0]
    
    if n != np.shape(A)[1]:
        raise Exception('Matrix is not square') 
    
    Q = np.eye(n)
    R = A

    for j in np.arange(0,n-1):
        # erzeuge Nullen in R in der j-ten Spalte unterhalb der Diagonalen:
        #a = (Q @ A)[j:,j:][:,0]
        a = np.copy((Q @ A)[j:,j:][:,0]).reshape(n-j,1)
        e = np.eye(n-j)[:,0].reshape(n-j,1)
        length_a = np.linalg.norm(a)
        if a[0] >= 0: sig = 1
        else: sig = -1
        v = a + sig * length_a * e          # vj := aj + sign(a1j) · |aj| · ej
        u = 1 / np.linalg.norm(v) * v       # uj := 1/|vj|*vj
        ut = u.T
        ua = u @ ut
        ub = 2 * ua
        x = np.eye(n)
        H = np.eye(n-j) - (2 * (u @ u.T))   # Hj := In − 2u1u1T bestimme die (n − j + 1) × (n − j + 1) Householder-Matrix Hj
        Qj = np.eye(n)
        Qj[j:,j:] = H                       # erweitere Hi durch einen Ii−1 Block links oben zur n × n Matrix Qi
        R = Qj @ R                          # R := Qj · R
        Q = Q @ Qj.T                        # Q := Q · QjT
        
    return(Q,R)


if __name__ == '__main__':
    # Beispiel aus Skript
    # A = np.array([[1, 2, -1], [4, -2, 6], [3, 1, 0]])
    # b = np.array([
    #     [9],
    #     [-4],
    #     [9]
    # ])

    # Beispiel aus Aufgabe 1
    A = np.array([
        [1, -2, 3],
        [-5, 4, 1],
        [2, -1, 3]
    ])
    b = np.array([
        [1],
        [9],
        [5]
    ])

    [Q,R]=Serie8_Aufg2(A)

    # Aufgabe 2b
    n = len(b) - 1
    QTb = Q.T @ b
    result = [0 for i in range(n+1)]
    row = n
    while row >= 0:
        value = QTb[row][0]
        column = n
        while column > row:
            value -= R[row][column] * result[column]
            column -= 1
        value = value / R[row,row]
        result[row] = value
        row -= 1

    print("\nQ:\n", Q)
    print("\nR:\n", R)
    print("\Result:\n", result)

    # Aufgabe 2c
    t1 = timeit.repeat("Serie8_Aufg2(A)", "from __main__ import Serie8_Aufg2, A", number=100)
    t2 = timeit.repeat("np.linalg.qr(A)", "from __main__ import np, A", number=100)
    avg_t1 = np.average(t1) / 100
    avg_t2 = np.average(t2) / 100

    print("Geschwindigkeit mit 3x3 Matrix:")
    print("Benötigte Zeit mit eigener Funktion:", avg_t1)
    print("Benötigte Zeit mit Numpy:", avg_t2)

    # Aufgabe 2d
    Test = np.random.rand(100,100)
    t1 = timeit.repeat("Serie8_Aufg2(Test)", "from __main__ import Serie8_Aufg2, Test", number=100)
    t2 = timeit.repeat("np.linalg.qr(Test)", "from __main__ import np, Test", number=100)
    avg_t1 = np.average(t1) / 100
    avg_t2 = np.average(t2) / 100

    print("Geschwindigkeit mit 100x100 Matrix:")
    print("Benötigte Zeit mit eigener Funktion:", avg_t1)
    print("Benötigte Zeit mit Numpy:", avg_t2)

    # Die von Numpy zur Verfügung gestellt Funktion ist wesentlich effizienter.