In our last installment, we discussed solutions to the secular equation. These solutions are the eigenvalues (and/or) singular values of matrices with a particular form. Since this post is otherwise light on technical content, I’ll dive into those matrix forms now. Continue reading

# DC SVD I: Just Can’t Let It Go

It’s been way, way, way too long since I’ve posted. I haven’t been slacking though, I’ve merely been busy. Really.

I decided to dive back into the SVD problem and look at an alternative to the QR based SVD computations. Namely, I’m going to give a breakdown of the *divide-and-conquer* approach to computing the SVD. Similar to the relationship between bidiagonal SVD and tridiagonal QR decompositions, there is a close relationship between dividing-and-conquering QR and SVD. I’m going to start from "the inside out" with the innermost task of DC (divide-and-conquer) process: solving a particular equation known as the secular equation (secular meaning "not heavenly" – i.e., earthly or planetly – check it out on wikipedia).

One note: I feel very guilty about still using Python 2. I had intented to transition this project to Python 3 over the course of the last year. Alas, my major training clients over the last year were using Python 2 and I didn’t really want to have a mixed development environment on my laptop. Well, at least there’s something to do if I make a book out of these posts!

Enough preamble, let’s get down to business. Continue reading

# SVD Computation Capstone

At long last, I’ve gotten to my original goal: writing a relatively easy to understand version of your plain, vanilla SVD computation. I’m going to spend just a few minutes glossing (literally) over the theory behind the code and then dive into some implementation that builds on the previous posts. Thanks for reading! Continue reading

# Householder Bidiagonalization

We’re going to make use of Householder reflections to turn a generic matrix into a bidiagonal matrix. I’ve laid a lot of foundation for these topics in the posts I just linked. Check them out! Onward. Continue reading

# Givens Rotations and the Case of the Blemished Bidiagonal Matrix

Last time, we looked at using Givens rotations to perform a QR factorization of a matrix. Today, we’re going to be a little more selective and use Givens rotations to walk values off the edge of a special class of matrices Continue reading

# Givens Rotations and QR

Today I want to talk about Givens rotations. Givens rotations are a generalization of the rotation matrix you might remember from high school trig class. Instead of rotating in the plane of a 2D matrix, we can rotated in any plane of a larger dimension matrix. We’ll use these rotations to selectively place zeros in a target matrix. Continue reading

# Householder Reflections and QR

My next two or three posts are going to deal with the QR and SVD techniques. QR can be done in several different ways. I’m going to work through two methods of computing QR that will shed some light on our SVD implementation. So, these two topics dovetail nicely. Onward! Continue reading

# Basic Cholesky Implementation

I spent a bunch of time talking about using lower level libraries (LAPACK directly and via LAPACKE or hand wrappers). My next set of posts is going to look at manually implementing these techniques. Mostly, I wanted to spend somem time with the numerical algorithms. I last visited them about 15 years ago in graduate school. There are also several performance lessons along the way (demonstating the limits of a "just compile it" for optimization). For starters, we’ll implement the Cholesky decomposition which is very easy to describe mathematically. Continue reading

# Testing a Cholesky Implementation

In working through the Cholesky, QR, and SVD methods for solving the normal regression equations, I got curious about the underlying implementations of each. I did QR and Cholesky in grad school (a great Numerical Analysis sequence). I never got around to SVD (I think we did some variation of tridiagonalization). Anyway, Cholesky is so simple that I think *testing Cholesky* generated more interesting points for me. There were two big issues and the "testing" post is (I think) longer than the "implementation" post.

Executing the tests brought up an interesting software engineering (code reuse) issue and creating good test inputs brought up an intersting mathematical issue. Read on! Continue reading

# Into the Weeds III – Hacking LAPACK Calls To Save Memory

In the last post, we tried to use lower-level LAPACK/LAPACKE calls to remove some of the copying overhead of our operations. Well, it turns out that it didn’t work out so well (see the table at the end of the post). The fundamental reason why is that LAPACKE took the easy way out, when they decided to wrap LAPACK functions. In particular, whenever there is a C-/Fortran-ordering issue (for a reminder), LAPACKE simply makes a copy. Continue reading