-1

I am trying to build function that performs QR decomposition, but upper triangular matrix R is not working, when I multiply it with Q the resulting matrix is different, I am not sure why it is not working properly.

import numpy as np
import math 

def GM(size, p,q,cos, sin):
    if (p >= size) or (p < 0) or (q >= size) or (q < 0):
            raise ValueError("Invalid first index i or j")
    matrixRotation = np.identity(size)
    matrixRotation[p,p]=cos
    matrixRotation[q,p]=-1 * sin
    matrixRotation[p,q]=sin
    matrixRotation[q,q]=cos
    return matrixRotation

def QRDecomposeWithGivens(matIn, *, traceOutput = True):
    m, n = matIn.shape
    matR = matIn.copy()
    matQ = np.identity(m)
    for col in range(n):
        for row in range(m-1, col, -1):
            x1 = matR[row-1, col]
            x2 = matR[row,col]
            d = math.sqrt(x1**2 + x2**2)
            cos = x1/d
            sin = x2/d
            rotationMatrix = GM(n, row-1, row,cos, sin)
            matR = rotationMatrix @ matR
            matQ = matQ @ np.transpose(rotationMatrix)
            if traceOutput:
                print("--- Next step ---")
                print("Rotation matrix for col {}, row {}; x{}: {}, x{}: {}".format(col, row, row-1, x1, row, x2))
                print(rotationMatrix)
                print("New R")
                print(matR)
                print("New Q")
                print(matQ)
                print("the result is")
                print (matQ*matR)
    return (matQ, matR)
    
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])
              



print(NEW_q*NeW_R)

QRDecomposeWithGivens(A)

I tried np.triu() function to build upper triangular matrix but it did not work

  • Looks okay to me, but you'd want to be using the `@` operator (for matrix multiplication, rather than element-by-element multiplication) to check your results rather than `*` – i.e. `print(matQ @ matR)` – motto Apr 30 '23 at 11:44
  • Welcome to Stack Overflow! Good places to get started here are the [tour] and [ask]. Could you show what output you're getting, and what output you expect? These are important parts of a [MRE]. – CrazyChucky Apr 30 '23 at 11:45
  • Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. – Community Apr 30 '23 at 14:57

1 Answers1

0

First, I would advise you to use the built-in function np.linalg.qr(A) for QR decomposition of matrix A, so you won't make a mistake. Second, your algorithm is correct. If you consider NEW_q, NeW_R = QRDecomposeWithGivens(A), then really A = NEW_q @ NeW_R. In your code, when you output the results, you don't use matrix multiplication, and so the resulting matrix is different from your expectations.

Mykola
  • 11
  • 2