Category: Math/Programming

Stuff about math and programming. Nerdy shit.

Posted on

Vectorization over C

The title is probably misleading, but this is a lesson I needed to talk about.

I wrote out some simple code for the quadrature over the reference triangle last time, which involves a double loop. To my chagrin, my immediate reaction to speeding up the code was to put it into Cython, and give it some type declaration.

This did speed up my integrals, but not as much as vectorization. By simply condensing one of the loops into a dot product, and using vector-function evaluation, I sped up my code a substantial amount, especially with higher order integration of “hard” functions.

def quadtriangle_vector(f, w00, w10, x00, x10):
    total = 0
    for i in range(len(w00)):
        total += w00[i] * np.dot(w10 / 2, f([(1 + x00[i]) * (1 - x10) / 2 - 1, x10]))
    return total

To see what I mean, consider the following function

from scipy.special import eval_jacobi as jac
def f(x):
    return jac(2, 1, 1, np.sin(x[0] - x[1]))
p = 20
x00, w00 = rj(p + 1, 0, 0)
x10, w10 = rj(p + 1, 1, 0)

The speedup I get is staggering.



Also, I tried to fully vectorize by removing the outer-loop. This actually slowed down the code a bit. Maybe I did it wrong? But for now, I’m decently happy with the speed.
Posted on

SPD Matrices and a False Inequality

This one is from my research, and it’s a doozy. Given two vectors x, y such that each element in x is less in absolute value than the corresponding element in y, show that for any SPD matrix A that x^TAx \le y^T A y.

After spending a good amount of timing trying to prove this, I realized that this is in general not true (in fact, the result I was suppose to be chasing would’ve been disproven if the above statement was true). As a counter example, consider the following counterexample from a Bernstein basis application:

Let x = [1, 0, -1/3], y = [1, -1, 1/3]. Let the matrix be

[2/7, 1/7, 2/35; 1/7, 6/35, 9/70 ; 2/35, 9/70, 6/35].

Then the quadratic forms will be 4/15 and 1/7 respectively.

Posted on

A Matrix Inequality

An exercise from Braess FEM book: for A, B symmetric, positive definite matrices, let A \le B. We want to show that the inverses satisfies a similar property B^{-1} \le A ^{-1}.

The book actually gave quite a lot of hints for this one.

\begin{aligned}x^TB^{-1}x &= x^T A^{-1/2} A^{1/2} B^{-1} x \\&\le \sqrt{x^T A^{-1} x} \sqrt{x^T B^{-1} A B^{-1} x}.\end{aligned}

.
Then from the hypothesis, A \le B \implies x^T(B - A)x \ge 0 \implies y^T(B^{-1} - B^{-1}AB^{-1})y \ge 0.
So we can plug this in to our equation above to find that x^TB^{-1}x \le \sqrt{x^T A^{-1} x}\sqrt{x^T B^{-1}x}.

Posted on

The (lack of a) Matrix

I think I finally understand why software packages like PETSc has an option for an operator when doing something like conjugate gradient. Why isn’t having a matrix good enough for everyone?

Well turns out that while all linear operators can be translated to a matrix, it may not be the best way to represent the operator. As an example, consider a basis transformation from Bernstein polynomials to Jacobi (or vice versa). It’s certainly possible to find and construct a matrix which does the operation, but it’s ugly.

On the other hand, it’s not that bad to write a code which utilizes the properties of the polynomials and convert it within in O(n^2) time. The key is that a Jacobi polynomial is a sum of Bernsteins, and Bernsteins can be degree raised or dropped at will.

This function will outperform the matrix in many sense. For one, there’s no need to construct a matrix, which will take n^2 operations in the first place. Next, matrix multiplication will take a n^3 operation, so if we optimize enough, we will always beat it. Finally, it’s really less painful to code, because each line of the function serves a visible purpose.

Anyways, I’m sold.

(I’ll eventually publish the code in the summer)

Posted on

Burgers Equation

For conservation laws in general, there’s the implicit solution of u = F(x - ut). It’s surprisingly useful for quick and dirty “exact” solutions to a lot of problems. In particular inviscid Burgers’ equation with initial values becomes trivial to code up.

For example, in Python it would be literally one line:

return optimize.fsolve(lambda u: u - self.u0(x - u*t), 1)
Posted on

Spectacular Books

Mathematicians suck at writing. This is part of the reasons why they went into math in the first place, because they suck at writing. Sucking at writing doesn’t mean that mathematicians don’t write books though; in fact there are tons of math books written by mathematicians.

The problem is most of them suck.

Half of them are for geniuses written by geniuses.

The other half are for geniuses written by borderline geniuses.

By far my two favorite books are William’s Probability with Martingales and Trefethen’s NLA book. Both British authors. Both written in a highly colloquial style.

It’s really too bad academics have this ego-stroking urge to write to the highest denominator rather than to, say grad students or undergrad.

God damn.

Posted on

Schur Complement and Minimal Energy Extension

(Note: this post is mainly for me to consolidate my thoughts)

In the framework of domain decomposition, consider creating the Schur complement which orthogonalizes interior nodes and the edges/vertex nodes. It turns out the norm of these functions which are orthogonal to the interior functions are minimal energy (i.e. L2 norm) extensions.

This can be seen in both a Hilbert space way or an optimization way. For the optimization way, write out the product for the mass matrix, and note that we can take a derivative to minimize one of the factors… now the Schur complement pops up naturally!

Posted on

Referencing Copies

A lot of times in numerical methods, I need to have a temporary variable as a time stepping tool or as an “incrementing” device without altering the original variable. A code snippet like

copy_u = u
for i in range(len(copy_u)):
    copy_u[i] = f(u)

But I gotta be more careful. There’s a deep difference between making a copy of a variable, and just creating a reference to a variable. Any change to copy_u might change u itself! Use np.copy or the copy module in Python!

Posted on

Bernstein Polynomials

This post won’t have too much math, but mainly musings.

After a four month hiatus from research, it seems my dive back into the world of unknown math is being squandered by my ineptitude. My documentation of the progress made, and the comments in my code seems to be greatly lacking. I’ve spent a good chunk of time trying to recall progress.

I think one of the cooler things Bernstein polynomials can prove is that polynomials are dense (uniformly too IRCC) in continuous functions on a bounded domain. Durrett has a proof in his book using a probabilistic view, which seems much cleaner than an algebraic/analysis proof such as here.

Posted on

Harder Every Decade

There’s a problem I solved awhile back which I was quite proud of:

Determine, with proof, the largest number that is the product of positive integers whose sum is 1976. (IMO 1976)

Note that this was at a high school olympiad level. Two days ago, my little brother brought me a state Mathcounts preparation packet. One of the problems was

Write 33 as the sum of two or more distinct prime numbers so that the product of these prime numbers is largest. What is the largest possible product?

While a proof is not needed, the move from an high school international level to a middle school state level in 40 years is pretty is pretty impressive.