Briefly, AI was able to generate from-scratch Python code for a function that computes the singular value decomposition (SVD) of a matrix — one of the most complicated functions I’ve come across in 40+ years of software engineering.
Where do I begin? I’ve looked at SVD for, literally decades. I learned early on that SVD is insanely complicated. It took dozens of mathematicians several decades of research to come up with a stable algorithm. The standard implementation (LAPACK) is over 10,000 lines of Fortran code.
I sat down one morning and fired up ChatGPT, and just for hoots, I asked it to generate from-scratch Python code for SVD. For the next several hours, ChatGPT gave me version after version that almost worked, but failed for various reasons. I told AI why the current version failed and it gave me a corrected version. Each new version got better and better.
The real point here is that each iteration from AI took only a few seconds to generate, instead of traditional development where each version would take hours or days.
In just a few hours, using AI, I was able to get a final working version of from-scratch SVD. Without AI, I estimate it would have taken me no less than 1,000 hours of full time effort.
Wow. Just Amazing.
The output of my demo:
Begin SVD demo TALL MATRIX A = [[ 3 1 1] [-1 3 1] [ 0 2 2] [ 4 4 -2]] U = [[-0.39318649 0.10749453 -0.80336030] [-0.21213925 -0.74140087 0.30624047] [-0.18505920 -0.60915594 -0.35957963] [-0.87530247 0.26018976 0.36267270]] s = [6.67747474 3.91438086 2.46758051] Vt = [[-0.66920959 -0.73394999 0.11608592] [ 0.53766957 -0.58611081 -0.60612338] [-0.51290346 0.34320770 -0.78685355]] reconstructed A [[ 3.00000000 1.00000000 1.00000000] [-1.00000000 3.00000000 1.00000000] [ 0.00000000 2.00000000 2.00000000] [ 4.00000000 4.00000000 -2.00000000]] ==================== np.linalg.svd results: U = [[-0.39318649 0.10749453 -0.80336030] [-0.21213925 -0.74140087 0.30624047] [-0.18505920 -0.60915594 -0.35957963] [-0.87530247 0.26018976 0.36267270]] s = [6.67747474 3.91438086 2.46758051] Vt = [[-0.66920959 -0.73394999 0.11608592] [ 0.53766957 -0.58611081 -0.60612338] [-0.51290346 0.34320770 -0.78685355]] reconstructed A [[ 3.00000000 1.00000000 1.00000000] [-1.00000000 3.00000000 1.00000000] [ 0.00000000 2.00000000 2.00000000] [ 4.00000000 4.00000000 -2.00000000]] WIDE MATRIX A = [[ 3 1 1 0 5] [-1 3 1 -2 4] [ 0 2 2 1 -3]] U = [[-0.72714903 0.19514593 -0.65815830] [-0.62191202 -0.59319511 0.51121913] [ 0.29065395 -0.78104906 -0.55270485]] s = [7.63921810 4.02386022 3.23278451] Vt = [[-0.20414852 -0.26332239 -0.10050154 0.20086846 -0.91571611] [ 0.29291099 -0.78196989 -0.48713106 0.10073440 0.23512158] [-0.76890186 -0.07111844 -0.38739015 -0.48724037 0.12750604]] reconstructed A [[ 3.00000000 1.00000000 1.00000000 0.00000000 5.00000000] [-1.00000000 3.00000000 1.00000000 -2.00000000 4.00000000] [ 0.00000000 2.00000000 2.00000000 1.00000000 -3.00000000]] ==================== np.linalg.svd results: U = [[-0.72714903 -0.19514593 -0.65815830] [-0.62191202 0.59319511 0.51121913] [ 0.29065395 0.78104906 -0.55270485]] s = [7.63921810 4.02386022 3.23278451] Vt = [[-0.20414852 -0.26332239 -0.10050154 0.20086846 -0.91571611] [-0.29291099 0.78196989 0.48713106 -0.10073440 -0.23512158] [-0.76890186 -0.07111844 -0.38739015 -0.48724037 0.12750604]] reconstructed A [[ 3.00000000 1.00000000 1.00000000 -0.00000000 5.00000000] [-1.00000000 3.00000000 1.00000000 -2.00000000 4.00000000] [-0.00000000 2.00000000 2.00000000 1.00000000 -3.00000000]] End demo
The code is from-scratch Python, but it uses a lot of NumPy calls such as np.diag(), and lots of matrix and array slicing. Maybe I’ll refactor the code to eliminate those calls by writing helper functions such as my_diag() but the code is very complex, and such a refactoring would take a big effort.
This was one of the most amazing technical explorations I have experienced in my lifetime.
Life is funny. At my age (mid-70s), the road ahead is shorter than the trail behind me. But I still wake up each day with a lot of joy in my heart, due mostly to family and friends but also unexpected things like the discovery in this blog post.
My last name, McCaffrey, is relatively rare so anyone with the name must be closely related. I came across an amusing headstone for one of my relatives. If you read the first letter of each sentence on the backside, you’ll see the real message. Apparently, John was quite the philanderer, and his wife met his mistress after John’s departure from this “mortal coil”. Ouch. Vulgar messages on headstones are not allowed in most cemeteries, so the wife and the mistress composed a message to John. Ouch.
Demo code. Replace “lt” (less than), “gt”, “lte”, “gte” with Boolean operator symbols (my blog editor often chokes on symbols).
# svd_from_scratch.py
import numpy as np
# -----------------------------------------------------------
def householder_qr(A):
"""
Reduced QR decomposition using Householder reflections
A: (m x n), m "gte" n
Returns Q: (m x n), R: (n x n)
"""
A = A.astype(float).copy()
m, n = A.shape
Q = np.eye(m)
R = A.copy()
for k in range(n):
x = R[k:, k]
normx = np.linalg.norm(x)
if normx == 0:
continue
v = x.copy()
sign = 1.0 if x[0] "gte" 0 else -1.0
v[0] += sign * normx
v /= np.linalg.norm(v)
# Apply reflector to R
R[k:, k:] -= 2 * np.outer(v, v @ R[k:, k:])
# Accumulate Q
Q[:, k:] -= 2 * np.outer(Q[:, k:] @ v, v)
return Q[:, :n], R[:n, :]
# -----------------------------------------------------------
def symmetric_qr_eig(A, tol=1e-12, max_iter=2000):
"""
Eigen-decomp of symmetric matrix using QR iteration
Returns eigenvalues and eigenvectors.
"""
A = A.astype(float).copy()
n = A.shape[0]
V = np.eye(n)
for _ in range(max_iter):
# QR decomposition (using our own QR)
Q, R = householder_qr(A)
A_new = R @ Q
V = V @ Q
# Check convergence (off-diagonal norm)
off = A_new - np.diag(np.diag(A_new))
if np.linalg.norm(off) "lt" tol:
A = A_new
break
A = A_new
eigvals = np.diag(A)
return eigvals, V
# -----------------------------------------------------------
def svd_qr_eig_numpy_only(A, tol=1e-12):
# SVD using only NumPy operations:
A = np.asarray(A, float)
m, n = A.shape
if m "gte" n:
# QR reduction
Q, R = householder_qr(A)
# Symmetric eigenproblem
RtR = R.T @ R
eigvals, V = symmetric_qr_eig(RtR, tol=tol)
# Sort descending
idx = np.argsort(eigvals)[::-1]
eigvals = eigvals[idx]
V = V[:, idx]
S = np.sqrt(np.clip(eigvals, 0, None))
nonzero = S "gt" tol
U = Q @ (R @ V[:, nonzero] / S[nonzero])
return U, S[nonzero], V[:, nonzero].T
else:
# Wide matrix: transpose trick
U, S, Vt = svd_qr_eig_numpy_only(A.T, tol)
return Vt.T, S, U.T
# ===========================================================
print("\nBegin SVD demo ")
print("\nTALL MATRIX ")
np.set_printoptions(precision = 8,
suppress=True, floatmode='fixed')
A = np.array([[3, 1, 1],
[-1, 3, 1],
[0, 2, 2],
[4, 4, -2]])
print("\nA = ");print(A)
# U, s, Vt = svd_qr_bidiagonal(A)
U, s, Vt = svd_qr_eig_numpy_only(A)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)
A_reconstructed = U @ np.diag(s) @ Vt
print("\nreconstructed A"); print(A_reconstructed)
print("\n====================")
print("\nnp.linalg.svd results: ")
U, s, Vt = np.linalg.svd(A, full_matrices=False)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)
A_reconstructed = U @ np.diag(s) @ Vt
print("\nreconstructed A"); print(A_reconstructed)
# ===========================================================
print("\nWIDE MATRIX ")
A = np.array([[3, 1, 1, 0, 5],
[-1, 3, 1, -2, 4],
[0, 2, 2, 1, -3]])
print("\nA = ");print(A)
# U, s, Vt = svd_qr_bidiagonal(A)
U, s, Vt = svd_qr_eig_numpy_only(A)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)
A_reconstructed = U @ np.diag(s) @ Vt
print("\nreconstructed A"); print(A_reconstructed)
print("\n====================")
print("\nnp.linalg.svd results: ")
U, s, Vt = np.linalg.svd(A, full_matrices=False)
print("\nU = "); print(U)
print("\ns = "); print(s)
print("\nVt = "); print(Vt)
A_reconstructed = U @ np.diag(s) @ Vt
print("\nreconstructed A"); print(A_reconstructed)
print("\nEnd demo ")
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2025 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2025 G2E Conference
2025 iSC West Conference
You must be logged in to post a comment.