Many areas of numerical analysis make use of kernel functions to describe interactions between pairs of vectors. In boundary integral methods for solving PDEs, Green’s functions describe how a delta forcing at one location influences the solution at another location. In the study of Gaussian Processes, covariance functions describe how pairs of random variables behave with respect to one another. When these kernel functions decay smoothly with distance between points (e.g. functions with a power of r in the denominator), they display approximate separability away from the origin. This manifests in their matrix discretizations as hierarchical rank structure, wherein off-diagonal blocks are numerically low rank. This talk will focus on techniques for factorizing these hierarchical matrices and updating the factorizations following updates to the underlying problem geometry. We’ll discuss a new technique for performing direct solves with systems where hierarchical matrices are rank-deficient subblocks. Finally we’ll see results from applying and parallelizing these methods in a shape optimization setting.