QR Analysis

output.txt

Matrix Rank: 9
QR is Solution? True
These a solution in the span of Q_2.
Diagonal of R (9); Rank of A (9)
Is non-zero diag of R equal to rank? (True)
Is x_1 a solution? True
Is x_1 == x_0? True
Is our new solution x_0? True
Is our new solution x_0? True
Is our new solution x_0? True
Is our new solution x_0? True

qr_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Setup\n",
    "import numpy as np\n",
    "import scipy,random\n",
    "import matplotlib.pyplot as plt\n",
    "np.set_printoptions(precision=2)\n",
    "\n",
    "#n>m (over-determined)\n",
    "m,n = (10,15)\n",
    "\n",
    "#Make exmaple ternary matrix\n",
    "A = (2*np.random.rand(m,n)-1).round()\n",
    "\n",
    "#Reduce Rank\n",
    "for i in range(n):\n",
    "    ids = set(range(n))-set([i]) #Remove current column\n",
    "    n_cols = 2#Number of columns in convex combination\n",
    "    col_ids = random.choices(list(ids),k=n_cols) #List of random col ids\n",
    "    cols = A[:,col_ids] #Get a list of columns\n",
    "    weights = np.random.rand(n_cols)\n",
    "    weights /= sum(weights) #Calculate convex combinations of weights\n",
    "    convex_combination = np.add.reduce([w*v for w,v in zip(weights.T,cols.T)])\n",
    "    A[:,i] = convex_combination\n",
    "\n",
    "print(\"Matrix Rank:\",np.linalg.matrix_rank(A))\n",
    "\n",
    "\n",
    "#Generate solution\n",
    "x = (2*np.random.rand( n ) - 1).round() #Ternary soultion\n",
    "y = A@x\n",
    "\n",
    "# QR pivoting\n",
    "rank = np.linalg.matrix_rank(A)\n",
    "Q,R,P= scipy.linalg.qr(A,pivoting=True)\n",
    "rank = r = sum(np.around(np.diag(R),12)!=0)\n",
    "P = np.eye(len(P))[:,P]\n",
    "Q_1 = Q[:,:rank]\n",
    "Q_2 = Q[:,:-rank]\n",
    "R_1 = R[:rank,:rank]\n",
    "R_2 = R[:,:-(rank-m)]\n",
    "# Get a single solution\n",
    "z_1 = scipy.linalg.solve(R_1,Q_1.T@y)\n",
    "padding = np.zeros((n-r)) #zero padding\n",
    "x_0 = P@np.hstack([z_1,padding]) #One solution\n",
    "print(\"QR is Solution?\",np.allclose(y,A@x_0)) \n",
    "#>>> QR is Solution? True\n",
    "\n",
    "if np.allclose(Q_2.T@y,0):\n",
    "    print(\"There exists a solution to this system\")\n",
    "else:\n",
    "    print(\"These a solution in the span of Q_2.\")\n",
    "\n",
    "print(f\"Diagonal of R ({r}); Rank of A ({np.linalg.matrix_rank(A)})\")\n",
    "print(f\"Is non-zero diag of R equal to rank? ({r==np.linalg.matrix_rank(A)})\")\n",
    "#>>> Is non-zero diag of R equal to rank? (True)\n",
    "\n",
    "\n",
    "# Trying to get soultion set\n",
    "z_2 = scipy.linalg.solve(R_1,Q_1.T@R_2)\n",
    "zeros = np.zeros((n-r,m-r)) #padding\n",
    "L = -P@np.vstack([z_2,zeros])\n",
    "\n",
    "x_1 = x_0 + L@(10*np.random.rand(L.shape[1])-5).round()\n",
    "while not np.allclose(y,A@x_1):\n",
    "    x_1 = x_0 + L@(10*np.random.rand(L.shape[1])-5).round()\n",
    "print(\"Is x_1 a solution?\", np.allclose(y,A@x_1))\n",
    "print(\"Is x_1 == x_0?\", np.allclose(x_0,x_1))\n",
    "\n",
    "#Show other solutions\n",
    "for i,col in enumerate(Q_2):\n",
    "    x_n = x_0 + L@col\n",
    "    if np.allclose(y,A@x_n):\n",
    "        print(\"Is our new solution x_0?\",np.allclose(x_n,x_0))\n",
    "\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

qr_analysis.py

#Setup
import numpy as np
import scipy,random
import matplotlib.pyplot as plt
np.set_printoptions(precision=2)

#n>m (over-determined)
m,n = (10,15)

#Make exmaple ternary matrix
A = (2*np.random.rand(m,n)-1).round()

#Reduce Rank
for i in range(n):
    ids = set(range(n))-set([i]) #Remove current column
    n_cols = 2#Number of columns in convex combination
    col_ids = random.choices(list(ids),k=n_cols) #List of random col ids
    cols = A[:,col_ids] #Get a list of columns
    weights = np.random.rand(n_cols)
    weights /= sum(weights) #Calculate convex combinations of weights
    convex_combination = np.add.reduce([w*v for w,v in zip(weights.T,cols.T)])
    A[:,i] = convex_combination

print("Matrix Rank:",np.linalg.matrix_rank(A))


#Generate solution
x = (2*np.random.rand( n ) - 1).round() #Ternary soultion
y = A@x

# QR pivoting
rank = np.linalg.matrix_rank(A)
Q,R,P= scipy.linalg.qr(A,pivoting=True)
rank = r = sum(np.around(np.diag(R),12)!=0)
P = np.eye(len(P))[:,P]
Q_1 = Q[:,:rank]
Q_2 = Q[:,:-rank]
R_1 = R[:rank,:rank]
R_2 = R[:,:-(rank-m)]
# Get a single solution
z_1 = scipy.linalg.solve(R_1,Q_1.T@y)
padding = np.zeros((n-r)) #zero padding
x_0 = P@np.hstack([z_1,padding]) #One solution
print("QR is Solution?",np.allclose(y,A@x_0)) 
#>>> QR is Solution? True

if np.allclose(Q_2.T@y,0):
    print("There exists a solution to this system")
else:
    print("These a solution in the span of Q_2.")

print(f"Diagonal of R ({r}); Rank of A ({np.linalg.matrix_rank(A)})")
print(f"Is non-zero diag of R equal to rank? ({r==np.linalg.matrix_rank(A)})")
#>>> Is non-zero diag of R equal to rank? (True)


# Trying to get soultion set
z_2 = scipy.linalg.solve(R_1,Q_1.T@R_2)
zeros = np.zeros((n-r,m-r)) #padding
L = -P@np.vstack([z_2,zeros])

x_1 = x_0 + L@(10*np.random.rand(L.shape[1])-5).round()
while not np.allclose(y,A@x_1):
    x_1 = x_0 + L@(10*np.random.rand(L.shape[1])-5).round()
print("Is x_1 a solution?", np.allclose(y,A@x_1))
print("Is x_1 == x_0?", np.allclose(x_0,x_1))

#Show other solutions
for i,col in enumerate(Q_2):
    x_n = x_0 + L@col
    if np.allclose(y,A@x_n):
        print("Is our new solution x_0?",np.allclose(x_n,x_0))