# ## MIT License

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#

# In[2]:
import numpy as np
import numpy.linalg as la
# In[3]:
A = np.random.randn(3, 3)
# Let's start from regular old (modified) Gram-Schmidt:
# In[4]:
Q = np.zeros(A.shape)
q = A[:, 0]
Q[:, 0] = q/la.norm(q)
# -----------
q = A[:, 1]
coeff = np.dot(Q[:, 0], q)
q = q - coeff*Q[:, 0]
Q[:, 1] = q/la.norm(q)
# -----------
q = A[:, 2]
coeff = np.dot(Q[:, 0], q)
q = q - coeff*Q[:, 0]
coeff = np.dot(Q[:, 1], q)
q = q - coeff*Q[:, 1]
Q[:, 2] = q/la.norm(q)
# In[5]:
Q.dot(Q.T)
# Now we want to keep track of what vector got added to what other vector, in the style of an elimination matrix.
#
# Let's call that matrix $R$.
#
# * Would it be $A=QR$ or $A=RQ$? Why?
# * Where are $R$'s nonzeros?
# In[6]:
R = np.zeros((A.shape[0], A.shape[0]))
# In[7]:
Q = np.zeros(A.shape)
q = A[:, 0]
Q[:, 0] = q/la.norm(q)
R[0,0] = la.norm(q)
# -----------
q = A[:, 1]
coeff = np.dot(Q[:, 0], q)
R[0,1] = coeff
q = q - coeff*Q[:, 0]
Q[:, 1] = q/la.norm(q)
R[1,1] = la.norm(q)
# -----------
q = A[:, 2]
coeff = np.dot(Q[:, 0], q)
R[0,2] = coeff
q = q - coeff*Q[:, 0]
coeff = np.dot(Q[:, 1], q)
R[1,2] = coeff
q = q- coeff*Q[:, 1]
Q[:, 2] = q/la.norm(q)
R[2,2] = la.norm(q)
# In[8]:
R
# In[9]:
la.norm(Q@R - A)
# This is called [QR factorization](https://en.wikipedia.org/wiki/QR_decomposition).
# ----------
# * When does it break?
# * Does it need something like pivoting?
# * Can we use it for something?