Latest Posts (20 found)

Monte-Carlo simulations

Monte Carlo simulations (or methods ) is the technique of applying randomness and the Law of large numbers to the solution of various scientific and engineering problems. One of its first documented uses was by Stanislaw Ulam and John von Neumann for nuclear weapon simulations after WWII [1] . In this post I want to provide examples of some simple uses of Monte Carlo simulations. We'll start with the classical example of calculating the value of \pi by throwing darts. Suppose we take a square board and inscribe a quarter of a circle into it. We then proceed to throw darts at the board and record whether each dart hits inside or outside the quarter circle. Having thrown many such darts, we calculate the ratio of the darts inside the circle to the total number thrown. Assuming our darts are distributed uniformly over the square, by the Law of large numbers this ratio should approach the ratio of areas of the quarter circle A_{circ} to the full square A_{square} . With a square side length of 1, we have: Here's a visualization: Change the number of samples (dart throws) and click "Run" to regenerate. The code is very simple - here's a slightly sanitized version: You'll notice that the estimate is relatively poor - even with 1000 samples - if you click "Run" several times, some numbers will be way off mark. While this method does estimate \pi , it's not a particularly good estimate. I find that running ~10 billion samples is necessary to estimate it to 4 digits after the decimal with reasonable reliability. In general, for independent trials like these, the typical sampling error decreases in proportion to 1/\sqrt{N} , where N is the number of trials. This means that halving the error requires four times as many samples. While the \pi estimation may seem whimsical, it's an example of an important class of problems to which Monte Carlo simulation is applied: numerical integration. Our simulation estimates the area under the quarter-circle curve, which is a definite integral. Many integrals are very difficult to solve analytically, and much research has been done in the area of numerical analysis to develop methods to calculate integrals. Monte Carlo methods are particularly useful for high-dimensional integrals, where other numerical methods can become prohibitively expensive. A common use of Monte Carlo methods is estimating complex combinatorial calculations. These often don't have analytical solutions, and enumerating all options is intractable due to the scale of the numbers involved. As an example, let's consider the game of SET . Each SET card has four attributes: And the goal is to find a "set" - three cards that are either all different or all the same for each attribute separately . As an example, here's a hand with a single set; see if you can find it [2] : And the next hand doesn't have any sets: Here's a question: given a freshly shuffled SET deck, what are the odds that the first 12 cards drawn will have no sets among them? This question is difficult to answer without using a computer. It's easy to calculate the number of ways to deal a 12-card hand from a deck of 81: But how many of these hands have no sets? Enumerating 70 trillion SET hands and checking each one can take quite a while, and there is no straightforward counting formula to answer this question. Some clever methods can be employed to leverage symmetries and other mathematical properties of SET to cut down this search space considerably. Donald Knuth himself worked on this problem and came up with a neat program ( setset-all on his programs page ) that found 2,284,535,476,080 such hands. Therefore, the answer to our question is: There's a 3.23% chance that a randomly drawn hand of 12 cards from a full deck of SET will have no set in it. Let's see how we can use a Monte Carlo simulation to answer this question with relatively small effort, without deep knowledge of the mathematical properties of SET that enable cutting down the search space Knuth-style. We can use the following pseudo-code: After running 10 million simulated draws, I got an answer of 0.0323, which matches the real answer very closely. The Monte Carlo approach lets us solve rather complicated problems in a very simple way. Suppose we want to answer the same question for a hand of 15 cards; this would blow up the search space considerably - there are about 100x more ways to select 15-card hands than there are to select 12-card hands. But for a Monte Carlo simulation, we adjust one small parameter and get a very reliable [3] answer (about 0.00037, in case you were wondering). One domain where Monte Carlo simulations are ubiquitous is projections for retirement portfolios. Suppose someone prepares to retire with a total sum of 1 million dollars in their portfolio; they'd like to be able to draw $30,000 a year from the portfolio for their living expenses. Would that work? There's a large number of factors to take into account when analyzing this question, but for simplicity let's focus on just two: portfolio return and inflation. We can run a naive estimate, assuming average values: suppose an average yearly portfolio return of 4%, and average yearly inflation of 2% [4] Let's denote our portfolio return as r=0.04 , and inflation as q=0.02 . Then the real return each year is: Starting with $1,000,000, at the end of the year we'll have $1,019,600 and then draw $30,000 for living expenses [5] , ending with $989,600. If we continue this way, the money runs out after ~55 years, which means that a person retiring at the age of 65 should be reasonably safe, right? But this is very simplistic ; assuming just average returns is risky, because they do a poor job of representing reality, and many factors have uncertainty. For example, the sequence of returns matters a lot; a bad year (-10%) followed by a great year (+18%) would still count as "4% on average" but produces significantly less money than two consecutive +4% years. Inflation is also unpredictable, and sometimes correlated with portfolio returns; there could be bad years of high inflation and low / volatile returns that can wreak havoc on a portfolio. As we add factors (variance in yearly draws, mixed portfolios of stocks, bonds, real estate, life expectancy, unexpected events, changing tax laws etc.), relying on a single average estimate becomes increasingly more fraught. This is why Monte Carlo simulations are very popular in this domain: by drawing from reasonable distributions based on historical data, a Monte Carlo simulation can easily run a million different scenarios and provide estimates: for example, what are the odds of money running out before death. Here's a useful chart from a simulation I ran: In the top chart: In the bottom chart: Looking a this simulation, under the current assumptions the plan sounds much riskier than the average assumption makes it appear. Assuming that a 65-y.o. person would plan for 25 years of retirement until death, the ~30% odds of not having sufficient funds for this duration of time are sobering. Perhaps a change in plans is needed (such as a more frugal lifestyle or securing additional funds in some way). Retirement projection is only one of may ways in which Monte Carlo simulations are used for financial and economical applications; given the high uncertainty of these domains, it's very difficult to plan using analytical calculations. Company sales projections, growth projections, stock offering prices and much more uses Monte Carlo simulations to arrive at estimates with reasonable error bars. All the code for the explorations in this post is available on GitHub . A quarter circle of radius 1 inside a unit square. Samples 200 10 1000 Run π ≈ — Change the number of samples (dart throws) and click "Run" to regenerate. The code is very simple - here's a slightly sanitized version: You'll notice that the estimate is relatively poor - even with 1000 samples - if you click "Run" several times, some numbers will be way off mark. While this method does estimate \pi , it's not a particularly good estimate. I find that running ~10 billion samples is necessary to estimate it to 4 digits after the decimal with reasonable reliability. In general, for independent trials like these, the typical sampling error decreases in proportion to 1/\sqrt{N} , where N is the number of trials. This means that halving the error requires four times as many samples. While the \pi estimation may seem whimsical, it's an example of an important class of problems to which Monte Carlo simulation is applied: numerical integration. Our simulation estimates the area under the quarter-circle curve, which is a definite integral. Many integrals are very difficult to solve analytically, and much research has been done in the area of numerical analysis to develop methods to calculate integrals. Monte Carlo methods are particularly useful for high-dimensional integrals, where other numerical methods can become prohibitively expensive. Combinatorial simulation - the game of SET A common use of Monte Carlo methods is estimating complex combinatorial calculations. These often don't have analytical solutions, and enumerating all options is intractable due to the scale of the numbers involved. As an example, let's consider the game of SET . Each SET card has four attributes: Number of shapes (1, 2 or 3) Color (Red, Green or Purple) Shape type (Oval, Diamond or Squiggle) Shading (Empty, Striped or Solid) The dashed line shows the constant assumptions mentioned before: what happens when yearly return is always 4% and inflation is always 2%. The shaded blue areas demonstrate the outcomes of 1,000,000 simulations where inflation and return numbers are drawn from reasonable normal distributions based on historical data. We see that in 25% of the cases, all money ran out by roughly 22 years. It's even easier to see how long the funds last; if we're interested in knowing, say, what are the odds that this plan will have enough money for 30 years - the chart shows it's about 60% (since in 40% of the simulations the funds were depleted at this point).

0 views

Summary of reading: July - September 2026

"Wuthering Heights" by Emily Brontë - good writing, but the protagonists are quite something. There's barely a likable character in the whole story - they are all either crazy, evil, stupid or a combination of these. It's also remarkable to consider how small the world of some people in those times was (pre-Victorian England); whole lives spent within a tiny area, interacting with a handful of other people. And that's for members of the nobility, not some bumpkin on a remote farm. "The Pickled City: The Story of New York pickles" by Paul Van Ravenstein and Monique Mulder - yes, a book about pickles, with lots of cool photos! That's what you get while randomly browsing in the local public library while your kids are picking out their books. All you wanted to know about the history of pickles, with particular emphasis on NYC in the early 20th century and the connection to Jewish immigration to that city. "Rebel Island: the incredible history of Taiwan" by Jonathan Clements - a detailed history of Taiwan, focusing mostly on the early years (pre-modern times). The writing is a bit academic and dry, but it's very comprehensive. I wish the book spent more time on the modern era, though. "The Old Ways: A Journey on Foot" by Robert Macfarlane - an ode to taking long walks (from hours to multi-day hikes). The writing is very beautiful, but much of it feels pointless; the author describes friends and obscure historical poets in a lot of detail which has little to do with the main thesis of the book. I'd say the first 1/4 of the good is very good, and the rest is unnecessary. "After Rain" by William Trevor - I'm on a quest to find additional authors that write evocative short stories that stay with you, and Trevor was recommended as someone to try since my current favorites are Steinbeck and Lahiri. "After Rain" is a collection of 12 stories, mostly about love and relationships - between partners but also parents and children; the settings are usually in 1980s Ireland or England. Very good stories overall, I'm glad I tried this author and should try something else by him. "Spoonhandle" by Ruth Moore - a novel about the lives of Maine natives living through changing times, seeing their picturesque coastal area become a magnet for "summer people" - rich vacationers from the big cities. There's no single protagonist in the book and the writing is very warm and human. I enjoyed it. "The Almighty Dollar: 500 Years of the World's Most Powerful Money" by Brendan Greeley - a history-oriented exploration of the importance of the dollar in world economy, from its beginnings as the Joachimsthaler in 16th century Bohemia and until modern times. This book has some interesting ideas but I found it to be quite academical and not easy to digest. "Behave: The Biology of Humans at Our Best and Worst" by Robert M. Sapolsky - the scope of this (very long) book is huge: explaining how human behavior works starting from basic biological principles. Our understanding of the actual biological underpinnings of behavior is so nascent, however, that the vast majority of the book focuses on the usual things - the psychological studies (including all the well known ones like prisoner's dilemma, the Stanford prison experiment, etc) and comparison with other animals. Having read several other books on this topic in the past, I haven't found much new in this one. "Humorous short stories" by Anton Chekhov (read in Russian) - a collection of short stories with a humorous angle. The stories are typically very short, just 2-3 pages each, but the character building is superb. I should try reading something longer by Chekhov. "Plainsong" by Kent Haruf - describes a year in the life of several loosely connected people in a remote farm town in eastern Colorado. A wholesome story, without a specific beginning or end - just the lives of ordinary folks. "Behind the Beautiful Forevers: Life, Death, and Hope in a Mumbai Undercity" by Katherine Boo - investigative journalism of the lives of several connected families in Annawadi - one of Mumbai's most miserable slums - roughly from 2007 to 2010. The only thing more shocking than the squalor in which these people live is the degree of corruption, fraud and exploitation surrounding them. A very poignant read. "The Dispossessed" by Ursula K. Le Guin - science fiction only in the sense of 1984; it's always amusing to examine how unambitious the futuristic vision of some authors is. The book was written in the early 1970s, and it presents life thousands of years in the future as one that's pretty much identical to what it was in the early 1970s, except for close-to-speed-of-light travel. But science is not the point of this book anyway; it's more of a political or philosophical fiction, examining an anarchist-communist society and comparing it to a capitalist, "propertarian" one. Overall, this book gave me very strong vibes of the famous Soviet kids story "Dunno on the moon" (Незнайка на Луне), which I recall reading as a kid. "The Lion Women of Tehran" by Marjan Kamali - a multi-generational novel focusing on the life of young women in Iran from the early 1950s to the early 1980s, in the context of the various authoritarian governments of those eras (from Mohammad Reza Pahlavi's rule to after the Islamic revolution of 1979). Nice book showing that a good story can be told without too many unnecessary twists and classical attention-grabbing dramas. "The Old Way: A story of the first people" by Elizabeth Marshall Thomas "Interpreter of Maladies" by Jhumpa Lahiri "Norwegian Wood" by Haruki Murakami "Cannery Row" by John Steinbeck

0 views

Rusty thoughts on "Parse, don't validate"

Like many programmers, I find Alexis King's Parse, don't validate article fascinating, because it gives a name to an idiom that seems familiar and important - one I've observed and used in the past without naming it explicitly. This post is a review of the "Parse, don't validate" pattern applied to the Rust programming language (the original post uses Haskell). I was particularly interested in finding educational examples of this pattern in the Rust standard library and other well-known projects. Without repeating the original article (please read it first!), here's the gist of it. Consider the venerable Vec ; its first method returns Option<&T> . Why? Because a vector is not guaranteed to have any elements in it, so what to do if first is invoked on an empty one? Returning an Option in this case is idiomatic in Rust [1] , with convenient syntax sugar for accepting the result of functions that return Option and deciding what to do next. So what's the issue? Imagine we have a function to read some configuration paths from an env var, while enforcing the invariant that the list can't be empty: So far, so good. Now let's take a typical usage of this function: Once get_configuration_directories returns a successful result, we are guaranteed that the vector isn't empty. And yet, if we want to get the first element of this vector, we have to use the first method that returns Option<&T> . We are therefore forced - again - to handle a potentially empty case (where the option is None ). As the original article states, this has a number of problems with code clarity, potential performance implications and a ticking time bomb if the invariant is ever changed in get_configuration_directories . The core issue is that Vec is fundamentally a type that can be empty; we can carry along a "This one can't be empty, pinky promise!" comment on all the relevant code, but it's not formally checked by anything. The solution is leveraging the type system to enforce a newly established invariant. We can use a separate type for "a vector that cannot be empty"; in fact, such types already exist in several Rust crates - for example nonempty : This type has no constructor that permits "no elements"; its new takes one element, and its first method returns &T without an Option : The rest of the crate deals with making NonEmpty behave as close as possible to a normal Vec , by implementing many useful traits, as well as conversions like: Let's see how our get_configuration_directories function would look if it returned a NonEmpty instead of a plain Vec : Note the use of NonEmpty::from_vec here - this is where the invariant is established. Now a successful result is NonEmpty , not just Vec . The client code looks like: There's no need to check if the returned value is empty again; this is enforced by the type system! This is where the parse vs. validate terminology of the original article comes from. When get_configuration_directories returned a Vec , it simply validated it. But when it returns a NonEmpty - the vector is transformed into another entity which carries additional meaning. If we treat the concept of parsing in the most generic sense - "transforming data from one format to another", this fits. To mention a less artificial example, the Rust rewrite of core POSIX utilities uses NonEmpty in several places [2] . For example, when constructing a shell pipeline: The command parser's code: A valid Pipeline is only returned if there are some commands in the parsed AST. Otherwise, it just returns None . Once this is done, the client code can use commands.first() without having to worry about the possibility of it returning None . A somewhat more interesting example can be found in the source code of rust-analyzer . This project has a type that represents an absolute filesystem path: Instead of carrying around a regular path, the absoluteness is recorded in the type once the initial parsing and validation is done: Subsequent code doesn't have to validate the the path is absolute. The type enforces it. Note also that AbsPathBuf wraps Utf8PathBuf , not PathBuf . Utf8PathBuf is itself a custom, "parsed" type refinement from the camino crate . Regular paths in the Rust standard library aren't guaranteed to be valid UTF-8, so they cannot be easily converted to a String (which has to be valid UTF-8 in Rust); camino::Utf8PathBuf establishes validity on construction, and can then be converted to a string with just: So we have an example of gradual parsing and type refinement here: Rust has a generic type called NonZero , to describe unsigned numeric quantities that are known to be non-zero. For example, thread::available_parallelism is defined as: If the call is successful, it returns a NonZero<usize> , which is like a normal usize with the restriction that it's not zero. Client code doesn't have to keep checking whether the parallelism is 0 - it's enshrined in the type system. Rust defines the division operator with NonZero<usize> in the denominator as an operation that "cannot panic". NonZero has an additional advantage: zero is an invalid value for the type, so Rust can use the zero bit pattern to represent None . Consequently, Option<NonZeroUsize> is guaranteed to have the same size and alignment as NonZeroUsize itself (and as usize ). This can avoid the extra storage that an Option<usize> would generally require. A common example of the "parse, don't validate" idiom appears in deserializing data from a JSON string. Rust's serde crate enables us to do the parsing, with validated decisions encoded into the type system, e.g.: And then later: There is a lot happening behind the scenes: We take code like this for granted these days, but it's still a great example of the pattern discussed in this post. Once the parser converted mode into the Mode enum, no further validation is required. In dynamic languages like Python and JS, the process is usually much more manual. Python's json.loads gives us a dictionary, and it's up to the user to validate its contents. Libraries like Pydantic permit an approach closer to Rust's, but they're not universally used. The types of all fields are enforced (e.g. "name" cannot be an array). The mode is validated to be one of the enum values of Mode . workers is validated to be a non-zero integer, because of the NonZeroUsize field type.

0 views

Notes on discrete-time Fourier series and transform

The following are my notes on discrete-time Fourier series (DTFS), as well as the discrete-time Fourier transform (DTFT). These topics serve as an important theoretical underpinning to the digital processing of signals by computers using the DFT (which will be covered in a future post). For discrete-time signals, we use the square bracket notation to denote samples: x[n] is the n-th sample of signal x . We’ll say that a discrete-time (just discrete from this point on, for brevity) signal is periodic with period N if: When talking about continuous-time Fourier series , we started with real trigonometric functions and later moved to complex exponentials. Here, we’ll just go ahead and start directly with (discrete) complex exponentials - converting between the two representations isn’t difficult and can be done as needed (Appendix C in this post demonstrates it). We’ll be dealing with the following family of signals [1] : These are discrete complex exponentials that are periodic with period N. w_0=\frac{2\pi}{N} is the angular frequency . As discussed in Appendix A, there are only N such distinct signals, since: We’ll want to consider the representation of an arbitrary N-periodic x[n] by a linear combination of \phi_k[n] : The notation k=\langle N\rangle means k runs over any sequence of N consecutive integers. Since there are only N distinct signals \phi_k , the order of summation doesn’t matter - as long as all of them are used. So it could go from 0 to N-1 , or from 1 to N , or from 2 to N+1 , and so on. All of these orders will enumerate all distinct \phi_k . Going back to our linear combination: This is the discrete-time Fourier series (DTFS) representation of x[n] and the coefficients a_k are the Fourier series coefficients. Note that there are no convergence issues as in the continuous case, because we’re dealing with a finite sum. To find the coefficients a_k , we’ll proceed in a way that’s somewhat similar to the continuous case. Multiply both sides of the equation above by e^{-ir(2\pi/N)n} and sum over N terms: Interchanging the order of summation on the right-hand side: By Appendix B, the inner sum on the right-hand side is equal to N when N\mid r-k and to 0 otherwise. Without loss of generality - since we iterate over N consecutive values - let’s just assume this happens in our sum when r=k (and therefore r-k=0 , which is an integer multiple of N ). So our equation becomes: To conclude, the Fourier decomposition of a periodic discrete signal x[n] into periodic complex exponentials is: Let’s revisit the triangular function t(x) from the post on Fourier series . Here, we’ll be using the following sampled version: The sampling is done with 12 samples per period of 4. The spacing between two samples is \Delta x=\frac{1}{3} . Another way to express it: Substituting this into the definition of t(x) , we get: We then do an odd extension [2] to the range n=[6..12) , and repeat it with period N=12 . Using considerations similar to the continuous time case , because our function is odd we just need a sine series here: The five terms correspond to the five pairs of nonzero-frequency indices; see Appendix C. And the coefficients are: We can run this sum over half the period, because the terms at and 12-n are equal (both t_d[n] and the sine change signs). Also, the terms at n=0 and n=6 vanish because t_d[0]=t_d[6]=0 : Looking at the values of t_d : We can similarly calculate all b_k until k=5 : The resulting Fourier series is therefore: You can use the following interactive plot to explore this approximation to the discrete triangle function t_d[n] : The dropdown box selects how many of the terms in the series shown above to plot. The orange crosses show where the series’ values at the integer indices are, and the orange line is interpolated to better visualize the sine waves being added. Note that when all the coefficients are used, the DTFS exactly reconstructs the input signal. In the post on the Fourier Transform we’ve seen how non-periodic signals can be represented in the frequency domain by taking the Fourier series and calculating them at the limit L\rightarrow\infty . Here, we’ll do something similar for discrete signals. Suppose we have a non-periodic discrete signal x[n] of finite duration (it’s zero outside a finite range of indices). Here’s a sample plot: The bottom half is x_N[n] - a periodic signal, one period ( N ) of which is equal to x[n] . We choose M large enough that x[n]=0 outside [-M,M] and set N=2M+1 . Since x_N is periodic, we can represent it with DTFS [3] : We’ll use the notation of angular frequency: Now the critical part; since x[n] is zero outside [-M,M] , extending the summation to all integers does not change its value: We’ll define the following continuous function with variable : Then at points w_k : The function X(w) is the discrete-time Fourier transform (DTFT) of x[n] . The inverse DTFT process reconstructs x[n] from X(w) , by substituting a_k back into the DTFS: By our definition of \Delta w , we have: Just like in the continuous case, we recognize the sum as a Riemann sum . Let M\rightarrow\infty . The spacing \Delta w tends to zero, and we can rewrite the sum as an integral: To conclude, the DTFT and inverse DTFT pair are: X(w) is a continuous, periodic function with period 2\pi , because: Therefore, it’s enough to run the integral on any interval of length 2\pi . The DTFS and DTFT are closely related. Starting with a finite-duration signal x[n] , we form a periodic signal by repeating a block of N samples containing it. Its DTFS coefficients are equally spaced, scaled samples of the DTFT: The relation is similar to the one between Fourier series and the Fourier transform in the continuous case. The DTFT also has useful properties similar to the continuous Fourier transform: linearity, time shifting, the convolution theorem etc, but we won’t spend time on them here. Discrete-time complex exponentials have the form e^{iw_0n} where w_0 is the angular frequency . The discrete nature of these functions leads to some interesting outcomes; for example, consider the exponential with angular frequency w_0+2\pi : Therefore, the exponential at w_0+2\pi is exactly the same as the exponential at w_0 . In considering discrete complex exponentials, we need only choose w_0 from a range of length 2\pi . Another interesting aspect of these functions concerns their periodicity. In order for e^{iw_0n} to be periodic with some positive integer period N , the following must hold: In other words: So w_0N itself must be an integral multiple of 2\pi . For some integer k : For a fixed positive integer N , let’s consider all complex exponentials for which N is a period: But earlier we’ve said that all complex exponentials with angular frequencies 2\pi apart are identical. This means that the set above consists of only N distinct exponentials. We can choose any starting k and get distinct signals with k , k+1 , k+2 and so on until k+N-1 . After that, they begin repeating: \phi_N[n]=\phi_0[n] and so on. All of this is very different from the continuous case, where t is a real number. The signals e^{iw_0t} and e^{i(w_0+2\pi)t} agree only at integer values of t , so they are distinct functions of t\in\mathbb{R} . Consequently, for a fixed positive period , the continuous-time exponentials e^{ik(2\pi/T)t} are distinct for every integer k , giving infinitely many distinct signals. Let’s consider the following discrete function, defined by a sum: We’ll want to consider two cases: (1) When k is an integer multiple of N : k=rN for some integer r (for example k=0 , k=N , k=-2N etc.) In this case, the sum becomes: Since both r and are integers, the exponent is an integer multiple of 2\pi i , so: (2) When k is not an integer multiple of N , let’s do a variable change p=n+1 : We can now use the finite sum formula for a geometric progression with the common ratio r=e^{i(2\pi/N)k} ; the sum is: Note that r\neq 1 because k is not an integer multiple of N . However: While these calculations demonstrate the sums from n=0 to N-1 , they work exactly the same for any N consecutive indices, because the summand is periodic in with period N . Therefore: Using the formula for DTFS developed earlier in the post: Recall that t_d[n] is odd over the full range n=[0..12) . Specifically, since the cosine is an even function, we have for any m : Therefore the cosine terms cancel out throughout the period; term 1 is the same as negative term 11, etc. We’re left with the sines: All coefficients are purely imaginary. For convenience, let’s write: Let’s turn back to the Fourier reconstruction formula: And pair up exponentials at k with ones at 12-k : Also, since t_d[n] are real, the conjugate of a_k is: So adding the paired exponentials with their coefficients and using Euler’s formula: We don’t include the terms at k=0 and k=6 , because a_0=a_6=0 (sine of an integral multiple of \pi ).

0 views

How big are factorials?

The other day, I found myself wondering how big 52! (52 factorial) is, and that led me to ponder how these could be estimated without a calculator or a computer. It turns out there’s some fairly interesting math behind being able to estimate the size (number of digits) of a factorial reasonably accurately. This post will start by stating how to do the estimate, and if you’re curious you can read on for the math background. Without further ado, the approximation is: As an example, let’s use my original question, by estimating this for 52! Well, 52 divided by e is... 20-ish? And \log_{10}(20) is about 1.3 [1] ; therefore our estimate comes out to: The real answer is 68, so this is very close! In estimates like this - when you’re dealing with enormous numbers - being off by a couple of digits usually isn't a big deal. The Gamma function for real n>0 is defined [2] as: This integral does not have an analytic expression in the general case, but it does have a very useful property that we can take advantage of. Let’s see what \Gamma(n+1) is: And now use integration by parts with: But notice that the last integral is just \Gamma(n) ; therefore, we’ve shown that: Let’s also calculate \Gamma(1) - it’s a special case that has an analytical solution: This helps establish an induction argument: In other words - the Gamma function is an interpolation of the factorial over all positive reals. Here’s a plot of the Gamma function over a small range; note that the y axis is log-scale because of the function’s fast growth: You may have encountered Stirling’s approximation before: It’s a great approximation that works reasonably well even for small . This section is a brief overview of how Stirling’s formula is derived from the Gamma function. We’ll start by massaging the integrand a bit: And making a change of variables x=ny , which means that dx=ndy : These steps make the integral amenable to applying Laplace’s method , which allows us to approximate definite integrals of the form: Where is a twice-differentiable function and some large number. By Laplace’s method, such integrals can be approximated by: Where is the global maximum of . Let’s see how to apply this method [3] to the latest equation we have for n! (renaming the dummy integration variable back to x ): In our case, f(x)=\ln x - x . It’s easy to show that this function is twice differentiable and has a global maximum at . Moreover: Substituting these into the proper places in Laplace’s approximation, we get: We can calculate the number of digits in n! by taking the base-10 logarithm of Stirling’s formula: Note that the first term is not multiplied by itself; therefore, as grows, it will become less and less noticeable. That said, it still adds a couple of digits - so you should take it into account if you want a more accurate approximation [4]

0 views

Concurrent Servers: Part 8 - Go

This is part 8 in a series of posts on writing concurrent network servers. In this part, we'll switch to Go and see how it tackles the challenges described earlier in the series. All posts in the series: This post assumes a basic familiarity with the Go programming language. As before, we'll start with a sequential server for the basic state machine protocol presented in part 1 . This is the main function: As in the previous parts, the server is "infinite"; it never stops serving new connections until it's explicitly killed. This is the function implementing the protocol for a single client; it takes a net.Conn value that represents a socket with a client connected on the other end: Rather than directly exposing OS threads, the Go runtime implements its own M:N scheduling of lightweight goroutines on top of OS threads. Using goroutines in Go is cheap - both in terms of syntax and developer effort, and in terms of system resources . Here's a version of our serial protocol server that serves clients concurrently by launching a goroutine for each client. The part of the code that's different from the previous sample is highlighted: The concurrent modification in this case is particularly simple because the server is infinite; there's no point waiting for these goroutines to finish (and hence no need for a sync.WaitGroup ). The parameters for server.ServeSerialProtocol are lexically captured from the enclosing scope and its return value is handled by the surrounding closure. Because goroutines are very cheap, this server is very unlikely to run out of resources due to launching too many goroutines; in fact, it will probably run out of something else - like file descriptors for sockets - first. However, sometimes it's still useful to limit the degree of concurrency - even in Go, and we'll discuss some approaches to do so in the following sections. Here are some scenarios in which it makes sense to limit the degree of concurrency in Go programs, even though goroutines are cheap to launch and operate: Let's switch to the primality testing server from part 4 for the rest of the post, because it represents a somewhat more realistic workload. As a reminder: the server receives numbers, simulates blocking by sleeping, and returns "prime" or "composite". The unbounded one-goroutine-per-client version looks almost identical to the previous code sample, except that the goroutine invocation calls another function: Where ServePrimeProtocol is [1] : The simplest way to limit concurrency in Go is by using a counting semaphore, implemented with a channel: The channel sem serves as a semaphore; note that it's a bounded channel with a maximal size. A token is acquired by sending to the channel, and released by receiving from the channel. When the channel is full, the send operation sem <- struct{}{} blocks until a token was removed by some other goroutine [2] . The type of the channel is struct{} which means "empty", or "no data". This is idiomatic in Go for channels that are used solely for their semantics, not to send/receive any actual data. Since launching goroutines is cheap and limiting concurrency is easy as shown above, the "worker pool" pattern is often unnecessary for scenarios like our server. Still, it has occasional uses (such as when workers have to maintain some non-trivial state across tasks) so it's worth discussing it here. Here's a variant of our primality testing server that uses a worker pool: A fixed number of worker goroutines is launched; these goroutines all receive "jobs" from the same channel. In the Accept loop, each client connection is sent to this channel as a new job and is picked up by the next available worker. As mentioned before, you would typically see a sync.WaitGroup somewhere to ensure clean shutdown of goroutines, but in our case it isn't necessary because we have a server that never exits. Do programmers have to resort to async / event-driven programming in Go? In my experience, almost never. Go was designed from the bottom up to be suitable for large-scale concurrency; goroutines are very cheap to create, have a tiny memory footprint and switching happens very quickly, all in user space. Measurements I ran back in 2018 have shown switching times of ~170 ns, as compared to 1-2 microseconds for threads on Linux. Moreover, Go already uses event-driven loops like epoll underneath for I/O. Goroutines that wait for I/O like sockets are effectively "parked" and consume no resources (beyond their small memory footprint); they are woken up by Go's runtime when their I/O descriptors are ready - this is very similar to how asynchronous programming works! That said, some people certainly do try to stretch their resources even more with direct asynchronous programming in Go when millions of streams are handled concurrently. All I'll say is that this is very rare, and an overwhelming majority of users never have to do this. In 2018, I wrote a post named Go hits the concurrency nail right on the head , and after several more years of active coding, I fully stand behind that statement. Go is extremely powerful and ergonomic for concurrent programs; while other environments go to great lengths to implement async-await style event loops in libraries, in Go it's already baked into the core language and runtime. You want event-driven I/O with very lightweight green threads that can also execute blocking tasks without worrying about the function coloring problem ? Go has you covered. All the code for this post is available on GitHub . Careful readers will note two issues with this code: (1) the protocol assumes the complete number is read from the socket in a single conn.Read call, and there's no framing - separation between distinct numbers; (2) the prime checking loop uses i*i which may overflow for large numbers. These issues are consistent across all the versions of the prime server in C, Python, JavaScript and Rust in earlier parts, because my focus was on the simplest possible code to demonstrate a point about concurrency. Part 1 - Introduction Part 2 - Threads Part 3 - Event-driven Part 4 - libuv Part 5 - Redis case study Part 6 - Callbacks, Promises and async/await Part 7 - Rust Part 8 - Go (this part) Tasks may be compute intensive, and the CPU capacity of any server is inherently limited. If too many concurrent goroutines compete for limited CPUs, they will all make very little progress. It may make more sense to have fewer tasks that complete in a reasonable time. Protecting potentially limited downstream resources, such as concurrent DB connections or other services. For example, if the server has to send requests to other services for each task, and these are rate-limited, concurrency will have to be carefully managed. Security reasons when work is dictated by clients; malicious clients can overload and crash a service that's too eager to serve, making it unavailable for legitimate clients.

0 views

Concurrent Servers: Part 7 - Rust

This is part 7 in a series of posts on writing concurrent network servers. In this part, we discuss how the challenges described in earlier parts are tackled in the Rust programming language. All posts in the series: Several years have passed since the previous parts were published. I've recently went over them to make sure the information presented is still relevant and all the code samples build and run using modern toolchains. I strongly recommend reviewing the previous parts before reading this one. This post assumes a basic familiarity with the Rust programming language. It will only explain Rust constructs when we encounter code that wouldn't appear in an introductory book or tutorial. The first few parts in the series focused on a socket server that implements a simple state machine protocol. See part 1 for a complete description of the protocol. Let's start by showing how this protocol is implemented in a basic sequential Rust server: With the function serve_connection defined as: As a reminder, this server version is sequential because it accepts clients one by one; the main loop blocks on serve_connection until it's done (the client closes the connection), and only then goes back to accept the next client. Clearly, handling clients one by one won't do. In part 2 , we've discussed approaches that use OS threads to handle clients concurrently. Let's start with the unbounded one-thread-per-client solution in Rust: The spawn method returns a Result<JoinHandle<T>> ; on success, we allow the handle to be dropped at the end of the loop iteration. In Rust, this detaches the thread; we don't actually wait for it to complete. This is reasonable for our code sample, because the loop is infinite ; it never terminates anyway. The potential for runaway threads is just one of the issues with the unbounded threads approach discussed in part 2. The solution is to use a fixed thread pool. Before diving into the code, a quick note on the design: the thread pool is a fixed set of threads that await "jobs" and handle them to completion. In our case a "job" is serve_connection for a specific client. There are many ways to implement a thread pool; for our use case, I went with a set of threads that all get a shared channel to which the main thread sends jobs. A worker thread picks up the next job from the channel, serves it to completion, and goes back to waiting for the next job. Here's how this looks in code: What is Receiver ? It's a type from the crossbeam_channel crate: Rust's builtin channels in std are mpsc - multi producer, single consumer, but what we need for our job queue is a channel that supports multiple consumers (the worker threads). While std does have mpmc , this is an experimental API only available in nightly versions at the time of writing. Therefore, I've opted to include the crossbeam_channel crate that provides well-tested mpmc channels for this sample [1] . And here's the main function: Note that our job channel is bounded - it has a fixed size. This helps naturally implement a backpressure mechanism - if too many clients connect, the following clients will have to wait - the main loop blocks on tx.send and won't accept additional clients on the socket until jobs are cleared from the channel. In parts 4, 5 and 6 of the series we've discussed event-driven , or asynchronous servers. Let's see how it's done in Rust. Specifically, part 6 presented a gradation from callbacks to promises to async/await mechanisms; Rust supports all of these and - as you'd expect - modern code is usually written with async/await while hiding all the details of promises (called futures in Rust) underneath. Without further ado, here's our simple state machine protocol in asynchronous Rust: Rust takes an interesting approach to async programming: it supports some of its fundamental building blocks (like futures and the async and await keywords) in the core language, but leaves the actual async engine implementation (the thing that implements the event loop) to external crates. By far the most popular crate for async programming in Rust is is tokio , so that's what we're using here. After reading the JS code in part 6, the Rust snippet above should appear fairly familiar, except perhaps the explicit tokio task "spawn". Instead of enqueuing a callback on the connection returned by listener.accept , the code spawns a tokio task, which can be seen as a green thread , and hence uses similar terminology [2] . These tasks must not issue blocking calls; therefore, they are supposed to use tokio's I/O utilities instead of the usual, blocking std utilities. In fact, we have to implement an async version of serve_connection to make this work: Note how similar this code is to serve_connection from earlier; the only real differences are the await calls on socket reads and writes [3] , and the types involved. For example, instead of a std::net::TcpStream used in the synchronous samples, here we're using tokio::net::TcpStream . Tokio has an underlying dependency called mio to handle non-blocking APIs for all kinds of I/O. It wraps OS-specific event loops like epoll to do so efficiently. While most of the series has been using a simple state machine server as the driving example, part 6 switched focus to a server for primality testing which simulates long compute tasks. Let's see how this is done in Rust with tokio: This code is very similar to the previous snippet conceptually; isprime is: Note that this sample demonstrates a job that can block (simulated with a sleep in this case). This can be problematic in an async context, as the tokio documentation explains . One potential solution would be to dispatch a blocking task to a separate thread pool and use tokio channels to communicate with it; this is similar to the approach we've taken in the thread pool sample above. Part 6 also included a version of this server that caches data on a local Redis instance; the goal was to demonstrate the complexity of event-driven code when additional layers of callbacks are added and how async/await can help mitigate that. Here's our Rust version of this server, using the redis crate (that has a tokio component enabled explicitly to support async calls): In conclusion, while Rust provides excellent support for async programming, it doesn't solve its inherent issues like function colors and the need for careful separation between blocking and non-blocking tasks. These issues are typically surmountable with some extra care, and async programming with Tokio in Rust is very popular due to its performance benefits. All the code for this post is available on GitHub . Part 1 - Introduction Part 2 - Threads Part 3 - Event-driven Part 4 - libuv Part 5 - Redis case study Part 6 - Callbacks, Promises and async/await Part 7 - Rust (this part) Because of the function color problem , the redis crate has a connection constructor specifically for async: get_multiplexed_async_connection . Here we have an example of shared state between tokio tasks - the Redis connection. Note that we don't require any particular synchronization because MultiplexedConnection is Clone ; cloning it to different tasks is safe - and in fact that's what we do for each new task. There's no magic here; if you look inside MultiplexedConnection , you'll see that it already has all the synchronization mechanisms implemented internally, as needed. Due to the magic of async/await, the code in serve_client is nice and linear. We simply await on the Redis call, and once it's back we continue with the rest of the handler. Since we're using an async Redis connection, in case waiting is required, control will be ceded to some other task that's not currently blocked on I/O.

0 views

Relative velocity and closing speed

In Physics simulations or game engines it’s sometimes useful to determine the speed with which two objects are approaching each other. This post will discuss the concept of closing speed , which is the normal component of the relative velocity of two objects. Suppose we have objects A and B [1] with velocity vectors \vec{V_A} and \vec{V_B} . The relative velocity of B w.r.t. A is: Put differently, it’s the velocity of B in A’s frame of reference. This relative velocity is a vector, and we can split it into orthogonal components. Obviously, the nature of such a split depends on the basis we want to use. We could look at the vector’s x an y components (we’ll be using - two dimensional space, but everything here applies to 3D as well), but for this post we’re interested in something slightly different: We draw a line connecting the two objects. The component of \vec{V}_{B|A} in the direction of this line is called the normal component of relative velocity, while the component perpendicular to this direction is called the tangential component . How do we find the component of a vector in the direction of a specific line? By using a vector projection! . We’ll represent the line by a vector, and find the projection of \vec{V}_{B|A} onto this vector. The positions of A and B can also be seen as vectors: \vec{P}_A and \vec{P}_B . The line connecting them can then be expressed as the vector \vec{P}_B-\vec{P}_A : All we need from this position difference vector is its direction, not its magnitude, however [2] . So we’ll use the unit vector of \vec{P}_B-\vec{P}_A , denoted as: Finally, to find the projection of \vec{V}_{B|A} onto \widehat{P} , we compute [3] : Where the multiplication operator between the vectors is the dot product . Note that the result of the dot product is a scalar; therefore, the quantity S_c is called the closing speed - it expresses the rate at which the relative distance of the two objects is changing. If it’s positive, the objects are drifting farther apart; if it’s negative, the objects are getting closer together. Therefore the term "closing speed" may be slightly confusing; alternatively, this has been called a "signed separation speed", or "normal relative speed" [4] . The signs in these calculations can be tricky to get right, so we have to be very careful. Let’s see a few examples that will help us make these computations more concrete. To build up some intuition and get some practice with the equations, we’ll review the following examples: Example I : We’ll start by computing the relative position unit vector \widehat{P} : Then, the closing speed is: Based on our convention, the negative sign of S_c means that the objects are approaching each other. Due to the simple nature of the example, this result is easy to verify, as it can be immediately guessed just by looking at the diagram. Example II : Here \widehat{P} is the same as in the previous example. The closing speed is: Same magnitude, but different sign from before, because the objects are moving farther apart. Example III This example is to demonstrate that we get consistent results even if B is to the left of A. Here the relative position unit vector is: And the closing speed: Which is the same as in example I, as expected. The direction of \widehat{P} flipped, but so did the direction of the relative velocity vector, so the result has the same sign. Example IV : Finally, an example showing more arbitrary positions and velocities. This example is a good opportunity to demonstrate something important about S_c : it’s time-dependent, because positions change with time. Here, -6.6 is the closing speed at the exact moment when A’s and B’s positions and velocities are as stated in the example. In the next time step, the position of A will be \langle2,5\rangle and the position of B will be \langle2,3\rangle , while their velocities remain the same. The S_c then will be quite different. This is a good segue to the next topic - which is a more physical view of closing speed. The computation shown so far represents a static view of the world; perhaps the right word to use is instantaneous . Given the positions and velocities at a given moment, what is the closing speed between the objects at that exact moment ? But there’s no reason to not generalize this using a more standard physical interpretation of velocity. First, let’s state the position vectors of A and B as a function of time: \vec{P}_A(t) and \vec{P}_B(t) . The relative position vector between the objects is also a function of time: Now we’ll define the scalar distance as the magnitude of this vector: We’re interested in \frac{dr(t)}{dt} - the change in this distance over time. Let’s start by breaking R(t) down to its constituents: By the chain rule [5] : Switching back to the vector representation: since \vec{R}(t)=\langle x(t),y(t)\rangle , the numerator of the fraction above is then a dot product between \vec{R}(t) and \vec{R}'(t) , we can write this as: But \hat{R}(t) is precisely \hat{P} from the earlier section, as a function of time. Moreover: Because velocity is the time derivative of position. Therefore, we end up with the same equation, just as a function of time: This formulation is more precise because it makes it very obvious that all the quantities we’re dealing with are time dependent.

0 views

Notes on the Fourier Transform

The Fourier series is a great tool for analyzing periodic functions. But what about functions that don’t repeat? We’ve seen that we can compute Fourier series for a non-periodic function defined on a finite interval, as long as we don’t care about its behavior beyond that interval. Let’s extend this idea to functions that never repeat; that is, non-periodic functions defined on the interval (-\infty,\infty) . To motivate the subject ahead, let’s look back at the example used in the earlier post about Fourier series : With an odd extension into [-2,0] . In that post, to make the Fourier series work, we assumed t(x) keeps repeating with a period 2L=4 on the entire x axis. Here, let’s face the reality that it does not - in fact - repeat, and observe how our Fourier series work out. Recall that the Fourier series approximating t(x) are the sine series (since it’s an odd function): The following visualization is interactive. By default, it shows t(x) (with its odd extension) and no Fourier series approximation. We’ll proceed by a series of steps and observe the outcome: Step 1 : set to some non-zero number; already at 3, the approximation is very good. The frequency spacing is \frac{\pi}{L} (this is the coefficient of x in the sines). Note that the Fourier series repeats every 2L , as expected. Step 2 : increase L to 6. This means our series are constructed assuming t(x) has a period of 12, not 4. Note how the Fourier series look now - they repeat every 12, and they don’t match t(x) as well as before. We can increase to a higher number to make the match better. As L grows, the spacing between adjacent frequencies decreases. Step 3 : increase L to 10. We no longer see the repetitions, so feel free to increase the values of x min and x max until you do. Note again that we need to add more and more coefficients to match t(x) better with this larger L , and the spacing adjacent frequencies grows smaller. Increasing L means our function repeats at larger and larger intervals. The logical conclusion of this progression is to ask - what happens if the function never repeats, meaning L\rightarrow\infty ? While not mathematically rigorous, the visual experiment here lets us make some conjectures: we’ll likely need an infinite number of coefficients for a good approximation, and moreover, the spacing between these coefficients will tend to zero. In other words, instead of a discrete set of coefficients, we’ll end up with a continuous line, or function . The function produced by this process is the Fourier transform of t(x) , and the next section shows its mathematical derivation. In these notes, we’ll be using the complex exponential formulation of Fourier series: We’re interested in a non-periodic defined on the interval (-\infty,\infty) . So we’ll be exploring the above equations for L\rightarrow\infty . First, let’s make a slight change of notation. Instead of writing formulae in terms of the period ( 2L ), we’ll be using the n-th harmonic angular frequency w_n : So we can slightly rewrite our series as: Using \Delta w as the difference between two consecutive frequencies: Using this notation, C_n is expressed as: So far there are no new insights here, just some new notation. Now we’re going to use it to facilitate the next step. Since L\rightarrow \infty , then \Delta w\rightarrow 0 . Let’s calculate the limit of the Fourier series representation of when \Delta w\rightarrow 0 : And substitute the latest C_n into this equation, changing its dummy integration variable from x to t to avoid confusion [1] Reordering slightly, and also replacing n\Delta w by w_n in the complex exponents: Looking at the limit with the sum carefully, this is a Riemann sum (see Appendix A)! w_n is the "sampled" version of , and \Delta w\rightarrow 0 . We can therefore replace it by an integral, changing w_n to and \Delta w to dw [2] : The inner integral is called the Fourier transform of and denoted [3] : And the full equation for is then the inverse Fourier transform: Let’s take our favorite odd triangular pulse example and calculate its Fourier transform. The function’s mathematical definition and plot are shown earlier in this post. Note that we’re not extending this function periodically - it’s zero beyond the range [-2,2] ; this is exactly why we need the Fourier transform here - as we’ve seen, Fourier series won’t do because the function they reconstruct eventually starts repeating. We’re looking to find: To calculate the integral, let’s decompose the complex exponent using Euler’s formula: Since our t(x) is odd, the first integral is zero . Also t(x)sin(wx) is even, so we can write: We’ve already calculated a very similar integral in the post on Fourier series , so let’s just skip to the result: The only remaining difficulty is its value at 0, which seems undefined at first (division by zero). However, note that as w\rightarrow 0 , the numerator also tends to 0, so we can use L’Hopital’s rule (twice!) to find that: This function is complex-valued; in fact, it’s purely imaginary. How do we visualize it? A common way to visualize complex-valued functions is by plotting their magnitude and phase separately. The magnitude of \hat{t}(w) is: Since \hat{t}(w) is purely imaginary, there are only two options for the phase: When the numerator is positive, we get a negative imaginary number with phase -\pi/2 , and when the numerator is negative, we get a positive imaginary number with phase \pi/2 . Finally, when \hat{t}(w)=0 (which happens at w=0 , by our earlier analysis, but also whenever is a whole multiple of \pi ), the phase is undefined. Here’s the magnitude and phase of \hat{t}(w) plotted against : It is common to talk about \hat{t}(w) as the frequency domain representation of t(x) . When the functions we’re working with have time as their domain (e.g. the x in t(x) represents time), which is often the case in the study of signals and systems, the Fourier transform can be seen as computing the frequency domain representation of the function. Here’s the Fourier transform formula again: It takes - the time domain representation of a function, and converts it to \hat{f}(w) - a frequency domain representation. For well-behaved functions, these two representations are dual - each one describes the function completely, just in a different way. To convert back from a frequency domain representation to the time domain, we use the inverse Fourier transform: While a time-domain plot ( t(x) ) shows how a signal changes over time, a frequency-domain plot ( \hat{t}(w) ) shows how the signal is distributed across all possible frequencies. Moreover, as we’ve seen, \hat{t}(w) is complex valued. Each frequency therefore has both a magnitude and a phase: the magnitude tells us how strongly that frequency contributes, while the phase tells us how that component is shifted. The frequency domain is extremely useful in signal analysis; for example, when designing filters. The Fourier transform also has a number of properties that are very useful in signal analysis and processing. But first, let’s discuss what a "well-behaved function" means for the purpose of applying Fourier transforms. The simplest existence condition for Fourier transforms is absolute integrability (also known as Lebesgue integrable): With this condition, \hat{f}(w) exists on the entire domain, is continuous and vanishes (tends to 0) as |w|\rightarrow\infty [4] . While this condition is sufficient, it’s not necessary; there are less well-behaved functions that also have Fourier transforms defined with some limitations. In these notes, we’re mostly interested in well-behaved functions that are used in real-world engineering, so we won’t discuss the other cases. Another assumption commonly made for real-world functions is that they vanish (tend to 0) as |x|\rightarrow\infty . While this is not a direct outcome of absolute integrability [5] , it’s a reasonable assumption in engineering. After all, real-world signals have finite energies. Intuitively, when we also assume is uniformly continuous , the assumption of vanishing at |x|\rightarrow\infty is a logical conclusion, because otherwise how can the total area for |f(x)| be finite? An important outcome of this discussion is that the Fourier transform is unsuitable for periodic functions. Functions that repeat at intervals are not absolute integrable . For periodic functions, we use Fourier series. The Fourier transform is a linear operator, because the integral is linear: So is the inverse Fourier transform; it’s similarly easy to show that: If we scale the domain of a function by a constant, its transform changes only slightly: Let’s do the variable substitution u=ax : This is the Fourier transform evaluated at \frac{w}{a} , so: There’s one small caveat here; when a is negative, the integral bounds should be flipped, causing a minus sign in front of the transform. So we can write: Which works for any a\ne 0 . This property is intuitive when thinking about signals: suppose a>0 , then f(ax) means the signal is compressed in the time domain by a factor a . The scaling property says that the frequency domain is expanded using the same factor; in other words, the higher frequencies become more prominent because we need sharper transitions to represent the compressed signal. Time shifting What happens to the Fourier transform if we time-shift the input signal by some constant: f(x-x_0) . By definition: Substituting u=x-x_0 , we get du=dx , so: Transform of a derivative An extremely useful property that’s often employed in the solution of partial differential equations; let’s calculate the Fourier transform of the derivative of : We’ll use integration by parts, where dv=f'(x) and u=e^{-i\cdot wx} . Therefore, v=f(x) and du=-iw\cdot e^{-i\cdot wx} : Recall the assumption made in the "Existence condition..." section about vanishing at infinities. So the first part of the equation above is zero, and we’re left with: Transform of convolution The convolution between two continuous functions and g(x) is defined as: Let’s calculate the Fourier transform of this function: This step of combining the integrals into a double integral, as well as the next step (changing the order of integration) is possible due to Fubini’s theorem and our assumption that and g(x) are Lebesgue integrable. Switch order of integration: Now, f(\xi) in the inner integral doesn’t depend on x , so we can pull it out: The inner integral is just the Fourier transform of a time-shifted g(x-\xi) , so we can write: And the remaining integral is the Fourier transform of , so: Convolution in the time domain translates to multiplication in the frequency domain! This result is so important in signal processing that it’s called the convolution theorem . Suppose we have some function and we want to know the area bounded between this function’s graph and the x axis in a certain interval [a,b] . One way to do this is to take a partition of the interval: And calculate the area under for every element of the partition. We can then approximate such sub-areas by rectangles, as follows: We’ll denote the area of each rectangle as f(x^*_i)\cdot\Delta x : There are many ways to choose which point of the interval [x_{i-1},x_i] to denote as x^*_i : left point ( x_{i-1} ), right point ( ), mid-point between the two (which is what our plot shows) or anything in between. The distinction doesn’t really matter for our purpose, as we will soon see. We can approximate the area under the curve of in the interval [a,b] with the Riemann sum , using a uniform partition: If is continuous on [a,b] , then as n\rightarrow \infty : This is known as the Riemann integral , or just the definite integral. The limit is why the exact choice of x^*_i doesn’t matter: as n\rightarrow\infty we have \Delta x\rightarrow 0 , and all points within [x_{i-1}, x_i] are equally good. The Fourier series is a great tool for analyzing periodic functions. But what about functions that don’t repeat? We’ve seen that we can compute Fourier series for a non-periodic function defined on a finite interval, as long as we don’t care about its behavior beyond that interval. Let’s extend this idea to functions that never repeat; that is, non-periodic functions defined on the interval (-\infty,\infty) . Visualizing Fourier series for non-repeating functions To motivate the subject ahead, let’s look back at the example used in the earlier post about Fourier series : \[t(x)= \begin{cases} x & 0 \leq x \leq 1 \\ 2-x & 1 < x \leq 2 \\ \end{cases}\] With an odd extension into [-2,0] . In that post, to make the Fourier series work, we assumed t(x) keeps repeating with a period 2L=4 on the entire x axis. Here, let’s face the reality that it does not - in fact - repeat, and observe how our Fourier series work out. Recall that the Fourier series approximating t(x) are the sine series (since it’s an odd function): \[t(x)=\frac{8}{\pi^2}\bigg[ sin\frac{\pi x}{2}-\frac{1}{3^2} sin\frac{3\pi x}{2}+\frac{1}{5^2}sin\frac{5\pi x}{2}-\cdots\bigg]\] The following visualization is interactive. By default, it shows t(x) (with its odd extension) and no Fourier series approximation. We’ll proceed by a series of steps and observe the outcome: n (terms in the Fourier series) L x min x max Step 1 : set to some non-zero number; already at 3, the approximation is very good. The frequency spacing is \frac{\pi}{L} (this is the coefficient of x in the sines). Note that the Fourier series repeats every 2L , as expected. Step 2 : increase L to 6. This means our series are constructed assuming t(x) has a period of 12, not 4. Note how the Fourier series look now - they repeat every 12, and they don’t match t(x) as well as before. We can increase to a higher number to make the match better. As L grows, the spacing between adjacent frequencies decreases. Step 3 : increase L to 10. We no longer see the repetitions, so feel free to increase the values of x min and x max until you do. Note again that we need to add more and more coefficients to match t(x) better with this larger L , and the spacing adjacent frequencies grows smaller. Increasing L means our function repeats at larger and larger intervals. The logical conclusion of this progression is to ask - what happens if the function never repeats, meaning L\rightarrow\infty ? While not mathematically rigorous, the visual experiment here lets us make some conjectures: we’ll likely need an infinite number of coefficients for a good approximation, and moreover, the spacing between these coefficients will tend to zero. In other words, instead of a discrete set of coefficients, we’ll end up with a continuous line, or function . The function produced by this process is the Fourier transform of t(x) , and the next section shows its mathematical derivation. Fourier series with L\rightarrow\infty leading to Fourier transform In these notes, we’ll be using the complex exponential formulation of Fourier series: \[f(x)=\sum_{n=-\infty}^{\infty}C_n\cdot e^{in\pi x/L}\] With: \[C_n=\frac{1}{2L}\int_{-L}^{L}f(x)e^{-in\pi x/L}dx\] We’re interested in a non-periodic defined on the interval (-\infty,\infty) . So we’ll be exploring the above equations for L\rightarrow\infty . First, let’s make a slight change of notation. Instead of writing formulae in terms of the period ( 2L ), we’ll be using the n-th harmonic angular frequency w_n : \[w_n=\frac{n\pi}{L}\] So we can slightly rewrite our series as: \[f(x)=\sum_{n=-\infty}^{\infty}C_n\cdot e^{i w_n x}=\sum_{n=-\infty}^{\infty}C_n\cdot e^{i\cdot n \Delta w x}\] Using \Delta w as the difference between two consecutive frequencies: \[\Delta w=w_n-w_{n-1}=\frac{n\pi}{L}-\frac{(n-1)\pi}{L}=\frac{\pi}{L}\] Using this notation, C_n is expressed as: \[C_n=\frac{\Delta w}{2\pi}\int_{-\pi/\Delta w}^{\pi/\Delta w}f(x)e^{-i\cdot n \Delta w x}dx\] So far there are no new insights here, just some new notation. Now we’re going to use it to facilitate the next step. Since L\rightarrow \infty , then \Delta w\rightarrow 0 . Let’s calculate the limit of the Fourier series representation of when \Delta w\rightarrow 0 : \[f(x)=\lim_{\Delta w\rightarrow 0}\sum_{n=-\infty}^{\infty}C_n\cdot e^{i\cdot n \Delta w x}\] And substitute the latest C_n into this equation, changing its dummy integration variable from x to t to avoid confusion [1] \[f(x)=\lim_{\Delta w\rightarrow 0}\sum_{n=-\infty}^{\infty}\left[\frac{\Delta w}{2\pi}\int_{-\pi/\Delta w}^{\pi/\Delta w}f(t)e^{-i\cdot n \Delta w t}dt\right]\cdot e^{i\cdot n \Delta w x}\] Reordering slightly, and also replacing n\Delta w by w_n in the complex exponents: \[f(x)=\frac{1}{2\pi}\lim_{\Delta w\rightarrow 0}\sum_{n=-\infty}^{\infty}\left[\int_{-\pi/\Delta w}^{\pi/\Delta w}f(t)e^{-i\cdot w_n t}dt\right]\cdot e^{i\cdot w_n x}\Delta w\] Looking at the limit with the sum carefully, this is a Riemann sum (see Appendix A)! w_n is the "sampled" version of , and \Delta w\rightarrow 0 . We can therefore replace it by an integral, changing w_n to and \Delta w to dw [2] : \[f(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\left[\int_{-\infty}^{\infty}f(t)e^{-i\cdot wt}dt\right]\cdot e^{i\cdot w x}dw\] The inner integral is called the Fourier transform of and denoted [3] : \[\boxed{\hat{f}(w)=\mathcal{F}\left[f(x)\right]=\int_{-\infty}^{\infty}f(x)e^{-i\cdot wx}dx}\] And the full equation for is then the inverse Fourier transform: \[\boxed{f(x)=\mathcal{F}^{-1}\left[\hat{f}(w)\right]=\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{f}(w)e^{i\cdot w x}dw}\] Example calculation of Fourier transform Let’s take our favorite odd triangular pulse example and calculate its Fourier transform. The function’s mathematical definition and plot are shown earlier in this post. Note that we’re not extending this function periodically - it’s zero beyond the range [-2,2] ; this is exactly why we need the Fourier transform here - as we’ve seen, Fourier series won’t do because the function they reconstruct eventually starts repeating. We’re looking to find: \[\hat{t}(w)=\int_{-\infty}^{\infty}t(x)e^{-iwx}dx\] To calculate the integral, let’s decompose the complex exponent using Euler’s formula: \[\hat{t}(w)=\int_{-\infty}^{\infty}t(x)cos(wx)dx-i\int_{-\infty}^{\infty}t(x)sin(wx)dx\] Since our t(x) is odd, the first integral is zero . Also t(x)sin(wx) is even, so we can write: \[\hat{t}(w)=-2i\int_{0}^{\infty}t(x)sin(wx)dx\] We’ve already calculated a very similar integral in the post on Fourier series , so let’s just skip to the result: \[\hat{t}(w)=-2i\cdot\frac{2\cdot sin(w)-sin(2w)}{w^2}\] The only remaining difficulty is its value at 0, which seems undefined at first (division by zero). However, note that as w\rightarrow 0 , the numerator also tends to 0, so we can use L’Hopital’s rule (twice!) to find that: \[\lim_{w\rightarrow 0} \hat{t}(w)=0\] Therefore: \[\hat{t}(w)= \begin{cases} -2i\cdot\frac{2\cdot sin(w)-sin(2w)}{w^2} & w\neq 0 \\ 0 & w=0 \\ \end{cases}\] This function is complex-valued; in fact, it’s purely imaginary. How do we visualize it? A common way to visualize complex-valued functions is by plotting their magnitude and phase separately. The magnitude of \hat{t}(w) is: \[|\hat{t}(w)|=\sqrt{\hat{t}(w)\cdot\hat{t}(w)^*}=2\left|\frac{2\cdot sin(w)-sin(2w)}{w^2} \right|\] Since \hat{t}(w) is purely imaginary, there are only two options for the phase: When the numerator is positive, we get a negative imaginary number with phase -\pi/2 , and when the numerator is negative, we get a positive imaginary number with phase \pi/2 . Finally, when \hat{t}(w)=0 (which happens at w=0 , by our earlier analysis, but also whenever is a whole multiple of \pi ), the phase is undefined. Here’s the magnitude and phase of \hat{t}(w) plotted against : It is common to talk about \hat{t}(w) as the frequency domain representation of t(x) . The frequency domain representation of functions When the functions we’re working with have time as their domain (e.g. the x in t(x) represents time), which is often the case in the study of signals and systems, the Fourier transform can be seen as computing the frequency domain representation of the function. Here’s the Fourier transform formula again: \[\hat{f}(w)=\mathcal{F}\left[f(x)\right]=\int_{-\infty}^{\infty}f(x)e^{-i\cdot wx}dx\] It takes - the time domain representation of a function, and converts it to \hat{f}(w) - a frequency domain representation. For well-behaved functions, these two representations are dual - each one describes the function completely, just in a different way. To convert back from a frequency domain representation to the time domain, we use the inverse Fourier transform: \[\mathcal{F}^{-1}\left[\hat{f}(w)\right]=\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{f}(w)e^{i\cdot w x}dw\] While a time-domain plot ( t(x) ) shows how a signal changes over time, a frequency-domain plot ( \hat{t}(w) ) shows how the signal is distributed across all possible frequencies. Moreover, as we’ve seen, \hat{t}(w) is complex valued. Each frequency therefore has both a magnitude and a phase: the magnitude tells us how strongly that frequency contributes, while the phase tells us how that component is shifted. The frequency domain is extremely useful in signal analysis; for example, when designing filters. The Fourier transform also has a number of properties that are very useful in signal analysis and processing. But first, let’s discuss what a "well-behaved function" means for the purpose of applying Fourier transforms. Existence condition for the Fourier transform The simplest existence condition for Fourier transforms is absolute integrability (also known as Lebesgue integrable): \[\int_{-\infty}^{\infty}|f(x)|dx<\infty\] With this condition, \hat{f}(w) exists on the entire domain, is continuous and vanishes (tends to 0) as |w|\rightarrow\infty [4] . While this condition is sufficient, it’s not necessary; there are less well-behaved functions that also have Fourier transforms defined with some limitations. In these notes, we’re mostly interested in well-behaved functions that are used in real-world engineering, so we won’t discuss the other cases. Another assumption commonly made for real-world functions is that they vanish (tend to 0) as |x|\rightarrow\infty . While this is not a direct outcome of absolute integrability [5] , it’s a reasonable assumption in engineering. After all, real-world signals have finite energies. Intuitively, when we also assume is uniformly continuous , the assumption of vanishing at |x|\rightarrow\infty is a logical conclusion, because otherwise how can the total area for |f(x)| be finite? An important outcome of this discussion is that the Fourier transform is unsuitable for periodic functions. Functions that repeat at intervals are not absolute integrable . For periodic functions, we use Fourier series. Some useful properties of Fourier transforms Linearity The Fourier transform is a linear operator, because the integral is linear: \[\begin{aligned} \mathcal{F}\left[\alpha f(x)+\beta g(x)\right]&=\int_{-\infty}^{\infty}\alpha f(x)e^{-i\cdot wx}dx+\int_{-\infty}^{\infty}\beta g(x)e^{-i\cdot wx}dx\\ &=\alpha\int_{-\infty}^{\infty}f(x)e^{-i\cdot wx}dx+\beta\int_{-\infty}^{\infty}g(x)e^{-i\cdot wx}dx\\ &=\alpha\mathcal{F}\left[f(x)\right]+\beta\mathcal{F}\left[g(x)\right] \end{aligned}\] So is the inverse Fourier transform; it’s similarly easy to show that: \[\mathcal{F}^{-1}\left[\alpha\hat{f}(w)+\beta\hat{g}(w)\right]= \alpha\mathcal{F}^{-1}\left[\hat{f}(w)\right]+\beta\mathcal{F}^{-1}\left[\hat{g}(w)\right]\] Scaling If we scale the domain of a function by a constant, its transform changes only slightly: \[\mathcal{F}\left[f(ax)\right]=\int_{-\infty}^{\infty}f(ax)e^{-i\cdot wx}dx\] Let’s do the variable substitution u=ax : \[\mathcal{F}\left[f(ax)\right]=\frac{1}{a}\int_{-\infty}^{\infty}f(u)e^{-i\cdot \frac{wu}{a}}du\] This is the Fourier transform evaluated at \frac{w}{a} , so: \[\mathcal{F}\left[f(ax)\right]=\frac{1}{a}\hat{f}\left(\frac{w}{a}\right)\] There’s one small caveat here; when a is negative, the integral bounds should be flipped, causing a minus sign in front of the transform. So we can write: \[\mathcal{F}\left[f(ax)\right]=\frac{1}{|a|}\hat{f}\left(\frac{w}{a}\right)\] Which works for any a\ne 0 . This property is intuitive when thinking about signals: suppose a>0 , then f(ax) means the signal is compressed in the time domain by a factor a . The scaling property says that the frequency domain is expanded using the same factor; in other words, the higher frequencies become more prominent because we need sharper transitions to represent the compressed signal. Time shifting What happens to the Fourier transform if we time-shift the input signal by some constant: f(x-x_0) . By definition: \[\mathcal{F}\left[f(x-x_0)\right]=\int_{-\infty}^{\infty}f(x-x_0)e^{-i\cdot wx}dx\] Substituting u=x-x_0 , we get du=dx , so: \[\begin{aligned} \mathcal{F}\left[f(x-x_0)\right]&=\int_{-\infty}^{\infty}f(u)e^{-i\cdot w(u+x_0)}du\\ &=e^{-iwx_0}\int_{-\infty}^{\infty}f(u)e^{-i\cdot wu}du\\ &=e^{-iwx_0}\mathcal{F}\left[f(x)\right] \end{aligned}\] Transform of a derivative An extremely useful property that’s often employed in the solution of partial differential equations; let’s calculate the Fourier transform of the derivative of : \[\mathcal{F}\left[f'(x)\right]=\int_{-\infty}^{\infty}f'(x)e^{-i\cdot wx}dx\] We’ll use integration by parts, where dv=f'(x) and u=e^{-i\cdot wx} . Therefore, v=f(x) and du=-iw\cdot e^{-i\cdot wx} : \[\mathcal{F}\left[f'(x)\right]=\left[f(x)e^{-i\cdot wx}\right]^{\infty}_{-\infty}-\int_{-\infty}^{\infty}f(x)(-iw\cdot e^{-i\cdot wx})dx\] Recall the assumption made in the "Existence condition..." section about vanishing at infinities. So the first part of the equation above is zero, and we’re left with: \[\begin{aligned} \mathcal{F}\left[f'(x)\right]&=-\int_{-\infty}^{\infty}f(x)(-iw\cdot e^{-i\cdot wx})dx\\ &=iw\int_{-\infty}^{\infty}f(x)e^{-i\cdot wx}dx\\ &=iw\cdot\mathcal{F}\left[f(x)\right] \end{aligned}\] Transform of convolution The convolution between two continuous functions and g(x) is defined as: \[(f\ast g)(x)=\int_{-\infty}^{\infty}f(\xi)g(x-\xi)d\xi\] Let’s calculate the Fourier transform of this function: \[\begin{aligned} \mathcal{F}\left[(f\ast g)(x)\right]&=\int_{-\infty}^{\infty}e^{-i\cdot wx}\left[\int_{-\infty}^{\infty}f(\xi)g(x-\xi)d\xi\right]dx\\ &=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-i\cdot wx}f(\xi)g(x-\xi)d\xi\ dx \end{aligned}\] This step of combining the integrals into a double integral, as well as the next step (changing the order of integration) is possible due to Fubini’s theorem and our assumption that and g(x) are Lebesgue integrable. Switch order of integration: \[\mathcal{F}\left[(f\ast g)(x)\right]=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-i\cdot wx}f(\xi)g(x-\xi)dx\ d\xi\] Now, f(\xi) in the inner integral doesn’t depend on x , so we can pull it out: \[\mathcal{F}\left[(f\ast g)(x)\right]=\int_{-\infty}^{\infty}f(\xi)\int_{-\infty}^{\infty}e^{-i\cdot wx}g(x-\xi)dx\ d\xi\] The inner integral is just the Fourier transform of a time-shifted g(x-\xi) , so we can write: \[\mathcal{F}\left[(f\ast g)(x)\right]=\int_{-\infty}^{\infty}f(\xi)e^{-i\cdot w\xi}\mathcal{F}\left[g(x)\right]d\xi=\mathcal{F}\left[g(x)\right]\int_{-\infty}^{\infty}e^{-i\cdot w\xi}f(\xi)d\xi\] And the remaining integral is the Fourier transform of , so: \[\mathcal{F}\left[(f\ast g)(x)\right]=\mathcal{F}\left[f\right]\cdot\mathcal{F}\left[g\right]\] Convolution in the time domain translates to multiplication in the frequency domain! This result is so important in signal processing that it’s called the convolution theorem . Appendix A: Riemann sum and the definite integral Suppose we have some function and we want to know the area bounded between this function’s graph and the x axis in a certain interval [a,b] . One way to do this is to take a partition of the interval: \[a=x_0<x_1<\cdots<x_{n-1}<x_n=b\] And calculate the area under for every element of the partition. We can then approximate such sub-areas by rectangles, as follows: We’ll denote the area of each rectangle as f(x^*_i)\cdot\Delta x : \Delta x=(b-a)/n is the width of one interval (assuming a uniform partition, but the math works just as well for non-uniform ones). x^*_i is some value in the interval [x_{i-1},x_i] .

0 views

Dot product: Component vs. Geometric definition

The goal of this post is to answer a simple question: why are the following two definitions of the vector dot product in Euclidean space [1] equivalent for vectors \vec{a} and \vec{b} : Here’s a graphical depiction of our vectors (focusing on for clarity, though this applies to any-dimensional vectors). It shows both the components of the vectors and the angle between them. The length of the arrow for \vec{a} is |\vec{a}| . We’ll show two proofs of the equivalence here, the geometric proof and the projection proof . The Appendix describes some properties of dot products that facilitate these proofs. We’ll be using this diagram of our vectors \vec{a} and \vec{b} , as well as the vector \vec{c}=\vec{a}-\vec{b} : Using the law of cosines [2] on the triangle formed by the three vectors: Since for any vector \vec{a} , we have \vec{a}\cdot\vec{a}=|\vec{a}|^2 (see Appendix), let’s rewrite this equation as: But \vec{c}=\vec{a}-\vec{b} and the dot product obeys the distributive property (see Appendix). Therefore: For this proof, we’ll assume the geometric definition is correct and will see how it leads to the component definition. We’ll begin by denoting vectors \vec{e}_1,\vec{e}_2\dots\vec{e}_n as the standard orthonormal basis for . For example, in 2D space, these basis vectors are \vec{e}_1=[1\ 0] and \vec{e}_2=[0\ 1] , shown in this diagram: If we take an arbitrary \vec{a}\in\mathbb{R}^n and calculate its dot product with a basis vector, we can use the geometric definition: where a_i is the component of \vec{a} in the direction of \vec{e}_i . The diagram makes it easy to see why this is true from basic trigonometry, but in the more general case this is just a vector projection . Now let’s represent vectors \vec{a} and \vec{b} as linear combinations of the basis vectors: And calculate the dot product \vec{a}\cdot\vec{b} , beginning by rewriting \vec{b} with its linear combination of basis vectors representation: Using the fact that the dot product distributes over linear combinations: But earlier we’ve shown that \vec{a}\cdot\vec{e}_i=a_i . Therefore: Which is the component definition \blacksquare . A generalization of dot products in is the inner product , which is an operation meeting some specific requirements, defined on a vector space. The inner product is denoted as \langle x,y\rangle:\mathbb{R}^n\times\mathbb{R}^n\to\mathbb{R} , and must satisfy the following requirements for all vectors x,y,z\in\mathbb{R}^n and scalars a,b\in\mathbb{R} : For , we define the inner product operation in its component formulation as: Let’s prove the requirements listed above for this operation; this is fairly straightforward, given the well-known properties of scalar multiplication and addition on : Linearity in the first argument: Positive-definiteness: Consider the components of vector x . Clearly, \forall i\quad x_i\cdot x_i=x_i^2\ge 0 . Since the vector x is not the zero vector, at least one of its components is nonzero, and for that component x_i\cdot x_i>0 . Therefore: Now that we’ve proved all the inner product requirements on our operation \langle x,y\rangle , we can say that is an inner product space with this operation. By meeting these requirements, it can be readily shown that our inner product operation has additional useful properties: The third property is particularly helpful, because it means the inner product is bilinear , and thus is distributive over addition. Note that these are shown for the component definition of dot product. It’s not too hard to prove distributivity for the geometric definition using the notion of projections and how they add up. The norm of a vector x in an inner product space is defined as |x|=\sqrt{\langle x,x\rangle} . Therefore, the square of the norm is |x|^2=\langle x,x\rangle . The norm is used to express the notion of magnitude , or length of a vector. If you think of a vector x\in\mathbb{R}^n in Cartesian coordinates, the definition of the norm is a generalization of the Pythagorean theorem. Component definition: \vec{a}\cdot\vec{b}=\sum_{i=1}^{n}a_i b_i Geometric definition: \vec{a}\cdot\vec{b}=|\vec{a}||\vec{b}|cos(\theta) , where |\vec{a}| is the magnitude of \vec{a} and is the angle between the vectors’ directions Symmetry: \langle x,y\rangle=\langle y,x\rangle Linearity in the first argument: \langle ax+by,z\rangle=a\langle x,z\rangle+b\langle y,z\rangle Positive-definiteness: if x\ne 0 then \langle x,x\rangle>0 \langle x,0\rangle=\langle 0,x\rangle=0 \langle x,x\rangle=0 if and only if x=0 \langle x,ay+bz\rangle=a\langle x,y\rangle+b\langle x,z\rangle \langle x+y,x+y\rangle=\langle x,x\rangle+2\langle x,y\rangle+\langle y,y\rangle

0 views

Summary of reading: April - June 2026

"The Nuremberg Trial" by John Tusa and Ann Tusa - a detailed, meticulously researched account of the Nuremberg Trials. There's not a whole lot of side questing in this book - it's all focused on the trials themselves. Interesting read overall, though somewhat dry and academic. "Things Become Other Things: A Walking Memoir" by Craig Mod - a kind of travelogue of the author walking across Japan's Kii peninsula, mixed with his childhood memories and impressions of life in Japan in general. It's a good book, though I thought I'd find more details about Japan here, whereas it's a much more introspective work about the author himself. "Social Justice Fallacies" by Thomas Sowell - the usual data-driven Sowell fare, using real historical data and statistical analysis to tackle some hot political issues like personal liberties, poverty data and affirmative action. "Focus: The ASML way" by Marc Hijink - I was inspired to read this book about ASML and its EUV technology after watching a fantastic Veritasium video on the topic. The book turned out to be a disappointment, however; I was interested in the technology behind ASML's machines, but the book is 98% focused on the human, political and organizational aspects of the company. If you're interested in the tech, the aforementioned video is a much better use of your time. "Every Living Thing: The Great and Deadly Race to Know All Life" by Jason Roberts - combined biography of Carl Linnaeus and the Comte de Buffon, who were groundbreaking naturalists in the 18th century working to categorize all living thiings. It's a history of the very early days of what was called "natural science", and later evolved into botany, biology, zoology, ecology and related disciplines. The title is hyperbolic, but the book itself is interesting and well written. "Junglekeeper: What It Takes to Change the World" by Paul Rosolie - the author is a conservationist in the Peruvian Amazonia rainforest. This book recounts his adventures on the path to establish Junglekeepers - his organization for preserving the forest and its wildlife. It's a nice read overall, though the writing is overly embellished and tiringly hyperbolic at times. "There Is No Place for Us: Working and Homeless in America" by Brian Goldstone - a poignant account of several families in Atlanta struggling with keeping access to suitable rental housing, circa 2020. It covers the danger zone of having just enough but absolutely no buffer, and how life emergencies affect families. Quite impressive investigative work to be able to produce a book like this; the stories are truly touching. The book does reasonably well to skirt around politics without too much preaching, which is appreciated. "Understanding Software Dynamics" by Richard Sites - from the ground up discussion of software performance analysis with sampling, tracing and understanding the different ways in which CPUs are unable to make progress. A good chunk of the book is dedicated to the author's KUtrace system. I wanted to like this book, but ultimately failed. I found it extremely verbose and tedious to follow, full of walls of text. "A Table for Two" by Amor Towles - a collection of short stories and a Novella. The stories all have some whimsical elements in them, and the writing is very good. That said, this book didn't quite recapture the magic of "A Gentleman in Moscow" for me. "Never Enough" by Jennifer Breheny Wallace - talks about the stress teens are under to excel academically and athletically to improve chances of admission to top universities, and what to do about it. Somewhat similar to "The Price of Privelege" and "Unequal Childhoods", though this one is more popular, in the sense that the author is a journalist, not a scientist, and the book is a collection of anecdata laced with official statistics, rather than experiences from the author's own research. "Thunder Below!" by Eugene B. Fluckey - the story of USS Barb, a submarine in the pacific during the latter part of WWII, written by its commander through 5 different deployments against Japan. Written in an engaging writing style, this book is very informative about how submarine warfare looked back then. I was somewhat shocked at the recklesness (suicidal courage?) of the commander though; he clearly was a very capable captain with a talented team, but surely there was lots of luck involved to be able to survive what he describes. That said, the book was also written 45 years after the events, so it's possible that there's some embellishment involved. "Breakneck: China's Quest to Engineer the Future" by Dan Wang - a very nice book about China's manufacturing superiority in the last few decades, as well as other aspects of its society like the one child policy and the COVID-19 lockdowns. The author contrasts China - "an engineer-driven society" with the USA - "a lawyer-driven society", discussing the effects on industrial capacity, culture and civil rights. Informative and well written. "Good People: A Novel" by Patmeena Sabit - an Afghan refugee family settles in Virginia in the early 2000s; this is a novel / mystery focused on their older children, their assimilation in the USA and what that lead to. Haunting book that will be difficult to get out of one's head, particularly for parents of teenage girls. "The Invention of China" by Bill Hayton - the author's thesis is that much of the national ethos of China - the unity of its peoples, language, territory - was invented about a century ago as part of a political agenda. The book is quite dry and academic, but this is key to its effectiveness, as it relies heavily on historical documents. It mentions a fantastic quote from Mao - "Make the past serve the present" - and I feel like this describes the book's main thesis very well. A fascinating example of the narrative of Orwell's 1948 in real life. "Advanced Hands-on Rust" by Herbert Wolverson - the idea of the book is to help one learn Rust through hands-on projects, by building games that use the Bevy framework. My conclusion is that Bevy is a particularly poor way to learn Rust because it's a massive, opaque and opinionated framework that bends your code to its will and conventions. Actually learning or practicing a language by building these simple games from scratch would be much better, IMO. So if your goal is to practice Bevy, sure, this book isn't too bad; but for Rust, stay away. Other than the Bevy issue, the book is poorly edited, with code samples out of sync with the accompanying code and diffuclt to follow to keep the project buildable. Also, since Bevy changes very quickly, you'll have to stick to the older version the book is using - otherwise things just won't compile. On the brighter side, the book does provide some coverage of Rust tooling that is useful - like benchmarking and creating well-behaved crates. But these topics in themselves are hardly worth a whole book. "Mathematics for Human Flourishing" by Francis Su - a math professor's attempt at defining the effects of doing mathematics on meaning in human life. I wanted to like this book, but unfortnately it's a bit too kumbaya for me. While I appreciate what the author was aiming at here, it just didn't click. "Benjamin Franklin: An American Life" by Walter Isaacson "A Gentleman in Moscow" by Amor Towles

0 views

Plugins case study: Pluggy

Recently I came upon Pluggy , a Python library for developing plugin systems. It was originally developed as part of the pytest project - known for its rich plugin ecosystem - and later extracted into a standalone library. You're supposed to reach out for Pluggy if you want to add a plugin system to your tool or library and want to use something proven rather than rolling your own. In this post I will share some notes on how Pluggy works, and will then review how it aligns with the fundamental concepts of plugin infrastructures . Pluggy is built around the concept of hooks : functions that host applications or tools (from here on, just "hosts") expose and plugins implement. A host exposes hooks by using a decorator returned from pluggy.HookspecMarker and a plugin implements this hook using a decorator returned from pluggy.HookimplMarker . Pluggy's documentation explains this fairly well; in this post, I'll show how to implement the htmlize tool with some plugins, introduced in the original article in my plugin series . As a reminder, htmlize is a toy tool that takes markup notation similar to reStructuredText, and converts it to to HTML. It supports plugins to handle custom "roles" like: As well as plugins that do arbitrary processing on the entire text. Out host defines two hooks: A hook is created by calling HookspecMarker with the project's name. This project name has to match between the host and its plugins. Pluggy is permissive about what hooks accept as parameters and what they return; for maximal flexibility and to stay true to the original htmlize example, our hooks return functions. To accompany this HookspecMarker , the host also defines a HookimplMarker with the same name: This is used by plugins to attach to hooks when they're loaded. The host's main function loads plugins at startup as follows: hookspecs is our Python module containing the hooks shown above. load_setuptools_entrypoints is Pluggy's helper for loading plugins that were pip -installed into the same environment and registered as setuptools entry points . It's a way to signal - in one's setup.py or pyproject.toml file - some metadata that projects can review at runtime. In our project, the plugins register themselves with this section in the pyproject.toml file: This says "for entry point htmlize , define a new entry named tt ". Pluggy's load_setuptools_entrypoints then uses importlib.metadata to access this information. Note that Pluggy doesn't require using this mechanism. Hosts can implement any plugin discovery method they want, and add plugins directly to their PluginManager with the register method. But this is the mechanism used for pytest and many other projects; it makes it very easy to automatically discover and register plugins that are installed with pip and equivalent tools. Once PluginManager loads the plugins, invoking them is straightforward; here's how htmlize invokes the contents hooks [1] : Generally, hook invocations return a list of all the hooks attached to by different plugins (a single host application can have multiple plugins installed and attaching to the same hook). When the host invokes the hook as shown above, the default order is LIFO, but plugins can affect this with hook options like tryfirst and trylast . Here's our entire narcissist plugin that's attaching to the contents hook: Some notes: Let's see how this case study of Pluggy measures against the Fundamental plugin concepts that were covered several times on this blog . It's important to remember that Pluggy is not a specific host application with a bespoke plugin system; rather, it's a reusable library for creating such plugin systems. Therefore, this is more of a meta case study. Generally, Pluggy leaves discovery logic to the user's discretion. Its PluginManager has a register method for adding plugins, and these can be discovered in any way the application chooses. That said, Pluggy comes with one discovery mechanism built in - through the entry points process of Python packaging, as shown above. This is hugely convenient for a large number of applications, as long as both the application and its plugins are installed via standard Python packaging tools (which is a very reasonable assumption in the Python ecosystem). In the entry point process, plugins register themselves by adding a [project.entry-points.<HOST-ID>] section in their pyproject.toml file. Otherwise - as in the previous section - users are free to devise their own registration schemes. This one is easy, since it's called hooks in Pluggy parlance as well! Pluggy's implementation of hooks is rather elegant, with function decorators available for plugins to set. We've seen an example of this above with @htmlize.hookimpl decorating htmlize_contents . Since Pluggy is designed for Python hosts and Python plugins, this one is fairly straightforward. The plugins typically assume the host project is already installed in the Python environment and its modules can be imported. In our example, hookimpl is imported from htmlize by the plugin to accomplish this. It also shows how host data is passed to the plugin - the post and db parameters. These are APIs exposed by the host for the plugins' use. In footnote 2 of my original fundamental concepts of plugin infrastructures post, I wrote [2] : I still believe my statement is true - plugin frameworks are very easy to create, and the functionality they provide is relatively small compared to their large surface area. In other words, this is a shallow API . That said, Pluggy does provide some nice functionality for the more advanced uses of plugins: Are these worthwhile for your project? It really depends on the project, and it's always worth keeping the tradeoff between dependencies and project effort in mind. The full code repository for this post is available here . It expects htmlize to be installed; as discussed previously, we rely on Pluggy's default install-based approach where both the host and plugins are installed into the same Python environment and can thus find each other. However, Pluggy supports any custom discovery method. It uses the hookimpl exported value shown earlier. It returns a function that acts on contents; this is the htmlize -specific contract (ABI, if you will) we've discussed before. Automatic entry point registration mechanism - if you need it Signature validation Consistent plugin result collection across multiple hook attachments in a single plugin and across many plugins Plugin ordering with firstresult , tryfirst , trylast , etc. Hook "wrappers" for some special use cases

0 views

Thoughts on starting new projects with LLM agents

A few months ago I wrote about using LLM agents to help restructuring one of my Python projects . It's worth beginning by saying that the rewrite has been successful by all reasonable measures; I've been able to continue maintaining that project since then without an issue. In this post, I want to discuss another project I've recently completed with significant help from agents: watgo . In this project many things are different; most notably, it's a from-scratch project rather than a rewrite, and it uses a different programming language (Go). This post describes my experience working on the project, and some lessons learned along the way. This is a new project, so it required extensive design. I began by iterating on the design with the agent, with a sketch of the API. For this purpose, I recommend using a Markdown file committed into the repository for future reference. After that, I started asking the agent to write CLs [1] in a logical order that made sense to me, keeping them small and reviewable (more on this in the next section). Sometimes it's not easy to have a small CL, and multiple rounds of revision may confuse the agent; in this case, I commit the CL and then go back and ask the agent to modify or refactor the code, as much as needed, with separate CLs. In the worst case, the whole sequence can be reverted if I feel we've taken the wrong direction (branches could also be helpful here for more complicated scenarios). This point is worth reiterating: sometimes a single CL is a huge step forward, but requires lots of review, cleanup and refactoring to be viable. I've had multiple instances where an agent produced several days of work in a single CL, but I then spent hours instructing it to clean up and refactor. Overall, it's still a productivity gain, just not as much as some pundits would like us to believe. Given the current state of agent capabilities, I think it's worth splitting projects into two categories: The watgo projects is a clear example of (2): I certainly intend to maintain this project in the long term, so I insist on code that I understand. With very few exceptions, no code gets in without full review and often multiple rounds of revisions. Even if the cost for writing code went down, maintaining a project is so much more than that. It's triaging and fixing bugs, it's thinking through what needs to be done rather than how to do it, it's keeping the code healthy over time, and so on. As Brian Kernighan said : Maybe at some point agents will become good enough that projects in category (2) can be implemented and maintained completely autonomously. Maybe. But we're certainly not there yet. My hunch is that getting there will require crossing the AGI line [2] , after which little in our world remains certain. If you're using an agent to send an actual PR and only review that , it's difficult to be disciplined enough to actually perform a thorough review. I find the following method to be more reliable: I use a CLI agent running locally in my repository, and ask it to update the code there. In parallel, I have a VSCode window open in the same project, where I can: Once I'm pleased with the change, I manually create a commit. As mentioned above, it's imperative to keep making progress in small chunks, with small enough CLs that a human can fully understand in a single review. It's very tempting to sprint ahead submitting thousands of lines of code every day, but this temptation has to be avoided. Coding with an agent is like speed-reading; yes, you're making more progress, but comprehension suffers the faster you go. Particularly for refactoring, agents still take the shortest route to destination. It's important to guide them to think about the "big picture" at all times, find all instances where X is better done as Y, not just a single place noticed during a review. This is why it's sometimes OK to have a CL submitted before you fully agree with everything, and go back to it later for several refactoring rounds. Source control works amazingly well when pair-coding with agents. It's a key point discussed in every "how to succeed with AI" article, but still critical enough to reiterate here: a solid testing strategy is absolutely crucial for success. Agents produce - by far - the best results when they have a solid test suite to test their code against. With the pycparser rewrite, I had a large existing test suite. For watgo , the very first thing I did was think through how to adapt the test suites of the WASM spec and of the wabt project for my needs. If your project doesn't have such tests to rely on, this should be your first order of business - finding one, or building one from scratch. Beware of self-reinforcing loops though; it's dangerous to trust agents for both the tests and the implementations tested against them. Go is a fantastic language for agents to write, because it's designed to be very readable by humans. The biggest strengths of Go are exactly what makes the experience of reviewing agent code so positive: Since most of the time spent by humans when using agents is reading rather than writing code, these effects compound and produce a great experience. Recall the discussion of how some languages are optimized for writability (Perl) while others are optimized for readability (Go)? Well, when working on a project with an agent we live in a world of 99% reading vs. 1% writing, so this really matters. I find this aspect really crucial in light of the earlier points made in this post - namely, keeping the human in the loop by understanding and reviewing all of the agent's design choices and code. If you're working on a subject that's completely new to you, I would strongly recommend against the approach described in this post. To really learn something, you have to work through it from scratch, yourself, reading, designing, writing the code. Agents don't change this basic fact; even before agents, if you wanted to learn X, copying it from Stack Overflow or some other project clearly wasn't the right way to go. Similarly, while agents can be used as a prop for learning, they cannot learn for you . As a corollary, junior engineers should exercise extreme caution when relying on LLMs. There's no replacement to hard-won experience and the sweat and tears of learning new, challenging topics. Learning is supposed to be hard; if it's too easy, you're probably not learning. For senior engineers, agents are a boon; it's a great tool to increase productivity, avoid the boring stuff, and get unstuck from procrastination; but only when used judiciously. Low importance / prototype / throw away projects where deep code understanding is unnecessary. These can be "vibe-coded" (submitting agent code without even reviewing it). High importance projects that I actually want to maintain; here, vibe-coding is ill advised and I insist on reviewing and guiding all code the agent writes before it's submitted (or shortly after, as discussed above). Review the agent's changes using VSCode's diff view Make my own tweaks and code changes if needed Go changes very infrequently, so you don't have to wonder "are we using the most modern / idiomatic approach" or "what the hell is this construct" as often as with other languages (looking at you, Python and TypeScript). There are relatively few ways to accomplish the same thing in Go, further lowering the mental burden. The standard library is rich and there's much less need to keep abreast of the package-everyone-uses du jour. In general, Go is designed for readability, with a mild-but-still-strong type system, uniform formatting, explicit error propagation and opinionated choices already made for you.

0 views

Notes on Fourier series

The trigonometric Fourier series is a beautiful mathematical theory that shows how to decompose a periodic function into an infinite sum of sinusoids. These are my notes on the subject, with some examples and the connection to linear algebra in Hilbert space. Let’s assume that is a well-behaved 2L -periodic [1] function and that we can find coefficients a_n and b_n such that: Then we say that the Fourier series on the right-hand side converges to . We’ll talk more about the assumptions mentioned above and convergence in the next section. Note that when n=0 , the sum becomes just ; therefore it’s customary to write the series starting with n=1 , with a separate constant component (which is the function's average over one period). To make computations nicer, this constant is typically called a_0 / 2 , so: Our goal is to find the coefficients a_n and b_n that satisfy this equation. We’ll do this in three steps. Step 1: Integrate both sides of the equation between -L and L [2] . Per Appendix A, all integrals within the sum are zero, so we’re left with: And thus we find : Step 2: Multiply both sides by cos\frac{m\pi x}{L} ( m is a positive integer constant) and integrate between -L and L . Looking at the right-hand side, the first integral is zero per Appendix A, and the last integral is zero per Appendix B. We’re left with: Per Appendix B, the integral on the right is zero for all n\neq m , and L for n=m . Therefore, we can write: Recall that m is an arbitrary integer, just like ; for consistency, we’ll replace m by and isolate a_n : Step 3: Hopefully it’s clear where this is going now; multiply both sides by sin\frac{m\pi x}{L} and integrate between -L and L . Using a very similar reasoning to step 2, we’ll end up with: We’ve just found a way to calculate all the coefficients of our Fourier series for : The previous section discusses Fourier series for a function that is well-behaved - but what does that mean? The full answer would lead us deep into analysis, which I’d like to avoid here. So I’ll keep it brief. We typically assume that is square integrable , which is denoted as L^2 . Moreover, we assume that the function is piecewise smooth : each segment of the function has continuous derivatives. A very simple example of a piecewise smooth function is f(x)=|x| . Another is the triangular wave function used in the example below. These conditions hold for pretty much any reasonable function we want to approximate using Fourier series, so they aren’t a serious burden. For a function that satisfies these conditions, it’s guaranteed to have a Fourier series that pointwise converges to it. This means that at every continuous point of , the Fourier series converges to it exactly; at every jump point, the Fourier series converges to the mid-point of the jump. Sometimes, additional properties of the function can help us simplify the Fourier series for it. If f_e(x) is an even function , then we know that: Because the function inside the integral is odd, and integrating an odd function over a symmetric interval results in 0. Therefore, the Fourier series for such f_e(x) is a cosine series : With coefficients and a_n given as before. Similarly if f_o(x) is an odd function, then its and a_n are 0, and its Fourier series is a sine series : So far we’ve been talking about 2L -periodic functions that can be faithfully represented by Fourier series. But what if we have a non-periodic function defined on a finite interval? E.g. suppose we have f(x)=x on the interval [0,L] . Can we approximate it with a Fourier series? Yes! First, we have to make a choice of how to extend the function to the negative interval [-L,0] . Then, we simply repeat the function every 2L - this is called a periodic extension . Note that the Fourier series calculation only cares about the range [-L,L] . The resulting series will approximate the generated periodic function in its entirety, and in particular will also converge to it in the [0,L] interval (except maybe the endpoints, depending on the mode of extension). There are several natural ways to extend a function defined on [0,L] into the interval [-L,0] [3] : Here’s an example of extending our sample function f(x)=x onto the full interval [-L,L] and then repeating it periodically every 2L : Note that the Fourier series for these extended functions will be different. However, they will all converge to in the interval [0,L] . Typically, even and odd extensions have the benefit of producing either cosine or sine series, correspondingly (as discussed in the previous section). We’ve seen that Fourier series work well for periodic functions and also non-periodic functions defined on a finite domain (because we can extend these periodically). But what about aperiodic functions defined on the entire real line? This is where we’ll have to leave Fourier series behind and move on to their generalization - the Fourier transform ; this will be a topic for a separate post. Let’s take the following triangular function t(x) [4] : t(x) is periodic with period 4. We can define it by starting with a formula on the interval [0,2] : Then making an odd extension into [-2,0] and repeating it periodically. Now we can go ahead to calculate its Fourier coefficients. Since this function is odd, we know that we’ll get a sine series , as a_n are going to be 0 for all . Let’s calculate b_n ; in our case L=2 (half the period). Since t(x) is odd and so is the sine, we’re integrating an even function over a symmetric interval. Therefore, we only have to integrate on the positive half of the range and multiply the result by two: Let’s set k=\frac{n\pi}{2} : And split up the integral for the different segments of t(x) : The first integral, by the method described in Appendix C: The second integral can also be split into two: The first of these is trivial to calculate; the second can once again use Appendix C. After some tedious but straightforward calculations [5] we’ll get: Adding I_1+I_2 , we get: Now let’s substitute k=\frac{n\pi}{2} back. This makes sin(2k) zero because the sine of an integer multiple of \pi is always zero: We have b_n , so the Fourier series for our t(x) is: Note that for even values of , sin \frac{n\pi}{2} is zero, so only the odd terms remain: Here’s an interactive chart showing how the series t(x) converges to our triangular function. You can set the number of terms in the Fourier series and see the effect (red line). Note that all even coefficients are zero so it will look the same for as for n-1 when is odd. We’ve written the Fourier series for as follows so far: We can rewrite this in a somewhat more compact form, using a single sinusoid with a configurable phase at each : Based on Appendix D, q_n and \theta_n can be computed as follows: When Fourier series are used in the context of signal processing, this formulation is easier to reason about because it represents the magnitude and phase shift of each harmonic of in the frequency domain [6] It should not come as a surprise that the Fourier series, being a combination of trigonometric functions, can also be represented with complex exponential functions. Specifically, we’ll show that our can be approximated as follows: Let’s calculate C_n . We proceed in a manner similar to before, by multiplying both sides of the equation by e^{-im\pi x/L} and taking an integral in the range [-L,L] : By Appendix A, the sum elements are all zero when n\neq m . When n=m , we get: Therefore, renaming m to (since it’s just an arbitrary integer constant): We’ve found an alternative formulation to Fourier series, using complex exponentials instead of trigonometric functions. While this was a direct derivation, another way to achieve the same result is to use the Euler Formula to derive: And substitute these into the original Fourier series formula. I’ll leave this as an exercise for the diligent reader; eventually, the result will be the same. Moreover, it’s possible to show a direct correspondence between a_n , b_n and C_n , for n>0 : Note that C_{-n}=C_n^* when both a_n and b_n are real (which is the case for a real-valued ). This helps explain why the complex formulation has negative frequencies in the sum; when the function is actually real, each negative frequency is paired up with a positive frequency and the result is real [7] : So, for a real function we only need to account for positive frequencies: We can take it further. C_n is a complex number, so let’s represent it in polar form as C_n=\frac{q_n}{2} e^{i\theta_n} (the factor of half will make sense soon). Then: And substituting back into the sum: This is precisely the compact formulation from the previous section! The most beautiful aspect of Fourier theory is that it doesn’t just happen to work by chance, and is deeply connected to linear algebra. Please read my post on Hilbert space before proceeding. The space of real-valued square integrable functions L^2 forms a Hilbert space, in which we can define the inner product (assuming real functions): We’ve demonstrated that the family of functions: Are all mutually orthogonal, because their pairwise inner products are zero! We’ve also shown that any function in L^2 can be represented as a weighted sum of these functions: So these functions form a basis for L^2 . When we think of these functions as vectors (in an infinite Hilbert space), much of what we did in this post starts feeling like "normal" linear algebra. For example, when we have a set of basis vectors and we want to know how to represent some vector in this basis, we usually find the coefficients by projecting it onto the basis. E.g. with a basis vector e_1 , the coefficient of : Similarly, when we calculate the coefficient b_n for some function , we project onto the basis vector sin\frac{n\pi x}{L} by calculating: From Appendix B, we know that the denominator is L , and we’ve just denoted: Which should look familiar! This is the core linear-algebra idea behind Fourier series: the functions 1 , cos\frac{n\pi x}{L} , and sin\frac{n\pi x}{L} play the role of orthogonal basis vectors, while the Fourier coefficients are coordinates of in this basis. The integral formulas for a_n and b_n are not mysterious tricks; they are projections, just like dot products with basis vectors in ordinary Euclidean space. Fourier series therefore let us decompose a function into independent orthogonal directions, much like decomposing a vector into its , , and z components. For any integer n\neq 0 and an arbitrary constant L, we have: Using these, we can calculate the integral of a complex exponential function for an integer n\neq 0 : We’ll start with the product of two sines, for any positive integers m and : Using the trigonometric identity for a product of sines, we can write: Now let’s focus on two different scenarios, m\neq n and m=n . If m\neq n , then each of the integrals constituting ss are 0 (see on Appendix A), so ss=0 . If m=n , then the second integral is still 0, but the first one isn’t: We can use exactly the same approach to show that: One more variant to cover: Since sine is an odd function and cosine is an even function, their product is an odd function. And the integral of an odd function over a symmetric interval is 0 (see this post for more details ). Let’s calculate the indefinite integral: For some constant k . We’ll use integration by parts: Here u=x , so du=dx . Also dv=sin(kx) , so v=-\frac{cos(kx)}{k} . Putting it together: Let’s take a general sinusoid with magnitude q , frequency and phase : We’re going to show that s(x) can be represented as a sum of a sine and a cosine with no phase. This is related to my earlier post on the sum of same-frequency sinusoids . Let’s start by expanding s(x) using a trigonometric identity: Now we’ll denote: a=q\cdot cos(\theta) and b=-q\cdot sin(\theta) , so: We have a and b in terms of q and , but what about the other way around? Let’s take the equations: Square both of them and add together: Now we’ll take the equations for b and a and divide one by the other: Where the atan2 function is careful to take into account the sign of both numerator and denominator. Also it’s worth mentioning that is determined up to additions of 2\pi . To conclude, for any q , and : With the aforementioned conversion formulas for a , b . The trigonometric Fourier series is a beautiful mathematical theory that shows how to decompose a periodic function into an infinite sum of sinusoids. These are my notes on the subject, with some examples and the connection to linear algebra in Hilbert space. Coefficients of Fourier series Let’s assume that is a well-behaved 2L -periodic [1] function and that we can find coefficients a_n and b_n such that: \[f(x)=\sum_{n=0}^{\infty}\left(a_n cos\frac{n\pi x}{L}+b_n sin\frac{n\pi x}{L}\right)\] Then we say that the Fourier series on the right-hand side converges to . We’ll talk more about the assumptions mentioned above and convergence in the next section. Note that when n=0 , the sum becomes just ; therefore it’s customary to write the series starting with n=1 , with a separate constant component (which is the function's average over one period). To make computations nicer, this constant is typically called a_0 / 2 , so: \[f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}\left(a_n cos\frac{n\pi x}{L}+b_n sin\frac{n\pi x}{L}\right)\] Our goal is to find the coefficients a_n and b_n that satisfy this equation. We’ll do this in three steps. Step 1: Integrate both sides of the equation between -L and L [2] . \[\int_{-L}^{L}f(x)dx=\int_{-L}^{L}\frac{a_0}{2}dx+\sum_{n=1}^{\infty}\bigg (\int_{-L}^{L}a_n cos\frac{n\pi x}{L}dx+\int_{-L}^{L}b_n sin\frac{n\pi x}{L}dx\bigg )\] Per Appendix A, all integrals within the sum are zero, so we’re left with: \[\int_{-L}^{L}f(x)dx=\int_{-L}^{L}\frac{a_0}{2}dx=\bigg[\frac{x\cdot a_0}{2}\bigg]_{-L}^{L}=a_0\cdot L\] And thus we find : \[a_0=\frac{1}{L}\int_{-L}^{L}f(x)dx\] Step 2: Multiply both sides by cos\frac{m\pi x}{L} ( m is a positive integer constant) and integrate between -L and L . \[\begin{aligned} \int_{-L}^{L}f(x)cos\frac{m\pi x}{L}dx&=\int_{-L}^{L}\frac{a_0}{2}cos\frac{m\pi x}{L}dx\\ &+\sum_{n=1}^{\infty}\bigg (\int_{-L}^{L}a_n cos\frac{n\pi x}{L}cos\frac{m\pi x}{L}dx+\int_{-L}^{L}b_n sin\frac{n\pi x}{L}cos\frac{m\pi x}{L}dx\bigg ) \end{aligned}\] Looking at the right-hand side, the first integral is zero per Appendix A, and the last integral is zero per Appendix B. We’re left with: \[\int_{-L}^{L}f(x)cos\frac{m\pi x}{L}dx=\sum_{n=1}^{\infty}\int_{-L}^{L}a_n cos\frac{n\pi x}{L}cos\frac{m\pi x}{L}dx\] Per Appendix B, the integral on the right is zero for all n\neq m , and L for n=m . Therefore, we can write: \[\int_{-L}^{L}f(x)cos\frac{m\pi x}{L}dx=a_m\cdot L\] Recall that m is an arbitrary integer, just like ; for consistency, we’ll replace m by and isolate a_n : \[a_n=\frac{1}{L}\int_{-L}^{L}f(x)cos\frac{n\pi x}{L}dx\] Step 3: Hopefully it’s clear where this is going now; multiply both sides by sin\frac{m\pi x}{L} and integrate between -L and L . Using a very similar reasoning to step 2, we’ll end up with: \[b_n=\frac{1}{L}\int_{-L}^{L}f(x)sin\frac{n\pi x}{L}dx\] We’ve just found a way to calculate all the coefficients of our Fourier series for : \[f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}\left(a_n cos\frac{n\pi x}{L}+b_n sin\frac{n\pi x}{L}\right)\] Where: \[\begin{aligned} a_0&=\frac{1}{L}\int_{-L}^{L}f(x)dx\\ a_n&=\frac{1}{L}\int_{-L}^{L}f(x)cos\frac{n\pi x}{L}dx\\ b_n&=\frac{1}{L}\int_{-L}^{L}f(x)sin\frac{n\pi x}{L}dx \end{aligned}\] Conditions on f and convergence of Fourier series The previous section discusses Fourier series for a function that is well-behaved - but what does that mean? The full answer would lead us deep into analysis, which I’d like to avoid here. So I’ll keep it brief. We typically assume that is square integrable , which is denoted as L^2 . Moreover, we assume that the function is piecewise smooth : each segment of the function has continuous derivatives. A very simple example of a piecewise smooth function is f(x)=|x| . Another is the triangular wave function used in the example below. These conditions hold for pretty much any reasonable function we want to approximate using Fourier series, so they aren’t a serious burden. For a function that satisfies these conditions, it’s guaranteed to have a Fourier series that pointwise converges to it. This means that at every continuous point of , the Fourier series converges to it exactly; at every jump point, the Fourier series converges to the mid-point of the jump. Cosine and Sine series Sometimes, additional properties of the function can help us simplify the Fourier series for it. If f_e(x) is an even function , then we know that: \[b_n=\frac{1}{L}\int_{-L}^{L}f(x)sin\frac{n\pi x}{L}dx=0\] Because the function inside the integral is odd, and integrating an odd function over a symmetric interval results in 0. Therefore, the Fourier series for such f_e(x) is a cosine series : \[f_e(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}a_n cos\frac{n\pi x}{L}\] With coefficients and a_n given as before. Similarly if f_o(x) is an odd function, then its and a_n are 0, and its Fourier series is a sine series : \[f_o(x)=\sum_{n=1}^{\infty}b_n sin\frac{n\pi x}{L}\] Fourier series for a non-periodic function defined on an interval So far we’ve been talking about 2L -periodic functions that can be faithfully represented by Fourier series. But what if we have a non-periodic function defined on a finite interval? E.g. suppose we have f(x)=x on the interval [0,L] . Can we approximate it with a Fourier series? Yes! First, we have to make a choice of how to extend the function to the negative interval [-L,0] . Then, we simply repeat the function every 2L - this is called a periodic extension . Note that the Fourier series calculation only cares about the range [-L,L] . The resulting series will approximate the generated periodic function in its entirety, and in particular will also converge to it in the [0,L] interval (except maybe the endpoints, depending on the mode of extension). There are several natural ways to extend a function defined on [0,L] into the interval [-L,0] [3] : Direct periodic repetition: we simply repeat every L : f(x+L)=f(x)\ \forall x . Even extension: f(|x|) Odd extension: when x\ge 0 and -f(-x) when x<0 .

0 views

Scaling, stretching and shifting sinusoids

This is a brief and simple [1] explanation of how to adjust the standard sinusoid sin(x) to change its amplitude, frequency and phase shift. More precisely, given the general function: We’ll see how adjusting the parameters , and affect the shape of s(x) . Each section below covers one of these aspects mathematically, and you can use the demo at the bottom to experiment with the topic visually. Scaling is conceptually the simplest change; we adjust to increase or decrease the amplitude (maximal height) of s(x) . Setting A=2 will make the value twice as large (in both the positive and negative direction) as the original function. Stretching changes the frequency of sin(x) , which is inverse proportional to its period. The baseline function sin(x) has a period of 2\pi , meaning it repeats every 2\pi . In other words, sin(x)=sin(x+2\pi) for any . If we set w=2 , we get sin(2x) . This function repeats itself twice as fast as sin(x) , because is multiplied by 2 before being fed into the sinusoid. If changes by \pi , the sinusoid’s input changes by 2\pi . Therefore, the period of sin(2x) is \pi , the period of sin(4x) is \frac{\pi}{2} and so on. [2] More generally, the period of sin(wx) is \frac{2\pi}{w} . Play with the demo below to see this in action, by changing and observing how the waveform changes. If we know the period p we want, we can easily calculate the that gives us this period: The final parameter we discuss is ; it’s called the phase of the sinusoid. In the baseline sin(x) , . The sinusoid is 0 at x=0 , achieves its positive peak at x=\frac{\pi}{2} , crosses 0 again at x=\pi , negative peak at x=\frac{3\pi}{2} and returns to its original position at x=2\pi where the repetition begins. By adding a non-zero , we don’t affect the sinusoid’s amplitude or frequency, but we do shift it right or left along the axis. For example, suppose we use the function sin(x+\theta) with \theta=\frac{\pi}{2} . Then when x=0 , we have sin(\frac{\pi}{2}) , so the sinusoid is already at its positive peak; at x=\frac{\pi}{2} , the sinusoid crosses 0 into the negatives, etc. Everything happens earlier (by exactly the value of \theta=\frac{\pi}{2} ) than in the baseline sinusoid. In other words, we’ve shifted the function left by \frac{\pi}{2} . Similarly, when is negative, everything happens later, and the function is shifted right . We’ve now gone over all the parameters for the function: Use the demo below to adjust these parameters and observe their effect on the sinusoid: controls the scaling factor (amplitude). is the frequency and controls the repetition period controls the phase - how much the sinusoid is shifted left or right

0 views

Thoughts on WebAssembly as a stack machine

This week the article Wasm is not quite a stack machine has been making the rounds and has caught my eye. The post claims that WASM is not a pure stack machine because it has locals and is missing some stack manipulation operations like dup and swap . While I don't necessarily disagree, IMHO it's a bit of a semantic discussion because - to the best of my knowledge - there is no formal definition of what is a stack machine. Wikipedia, for example, says: WASM certainly fits this definition; the primary interaction is through the stack, though WASM is augmented with an infinite register file (locals). The more purist stack machines like Forth are only limited to the stack and a memory (pointers into which are managed on the stack); WASM has these too, plus the registers. Speaking of Forth, the mention of dup reminded me of my own impressions of programming in that language, documented in my post about implementing Forth in Go and C . There, I highlighted the following essential library function for Forth; it adds an addend to a value stored in memory. And lamented how difficult it is to understand such code without the detailed stack view in comments alongside it. I find it much simpler to reason about this WASM code: You may say this is cheating because folded WASM instructions help readability and they're just syntactic sugar; OK, here's the linear code: It's still very readable, because - while the stack is used for all the calculations and actual commands - some of the data lives in named "registers" instead of on the stack. So we don't need all those tuck-swap contortions to get things into the right order. One might worry about the duplicated local.get $addr ; wouldn't a real dup be better? Well, not in terms of readability, as we've already discussed. How about performance? Since the stack VM is just an abstraction and the underlying CPUs executing this code are register machines anyway, the answer is no - it doesn't matter at all. Modern compiler engineers were forged in the fires of C and its descendants; arbitrary control flow, arbitrary register and memory access, anything goes. Compilers are quite sophisticated. Let's see how wasmtime compiles our add_to_byte to native code (using wasmtime explore with its default opt-level=2 ); comments are added by me: This is pretty much the code we'd expect to be emitted for the C statement mem[addr] += addend , or if we were writing x86-64 assembly by hand. The compiler had no difficulty figuring out that two consecutive loads from the same WASM local produce the same value and do not - in fact - have to be duplicated. The WASM model makes it rather easy, because you can't alias locals; as long as there are no intervening writes into the same local, multiple reads are known to produce the same value (redundant load elimination).

0 views

Debugging WASM in Chrome DevTools

When I was working on the WASM backend for my Scheme compiler , I ran into several tricky situations with debugging generated WASM code. It turned out that Chrome has a very capable WASM debugger in its DevTools, so in this brief post I want to share how it can be used. I'll be using an example from my wasm-wat-samples project for this post. In fact, everything is already in place in the gc-print-scheme-pairs sample. This sample shows how to construct Scheme-like s-exprs in WASM using gc references and print them out recursively. The sample supports nested pairs of integers, booleans and symbols. To see this in action, we have to first compile the WAT file to WASM, for example using watgo : The browser-loader.html file in that directory already expects to load gc-print-scheme-pairs.wasm . But we can't just open it directly from the file-system; since it loads WASM, this file needs to be served with a local HTTP server. I personally use static-server for this, but you can use anything else - like Python's built-in http.server : Now it can be opened in the browser by following the printed link and selecting the browser-loader.html file. Open the Chrome DevTools, and in Sources , open the Page view on the left. It should have one entry under wasm , which will show the decompiled WAT code for our module. Note: this code is disassembled from the binary WASM, so it will lose some WAT syntactic sugar (like folded instructions): You can set a breakpoint by clicking on the address column to the left of the code, and then refresh the page. The DevTools debugger will run the program again and stop at the breakpoint: Here you can step over, into, see local values and call stack, etc - a real debugger! The most important use case for me while developing the compiler was debugging unexpected exceptions (coming from instructions like ref.cast ). Notice the checkboxes saying "Pause on ... exceptions" on the right-hand side of the previous screenshot. With these selected, the DevTools debugger will automatically stop on an exception and show where it is coming from. Let's modify the gc-print-scheme-pairs.wat sample to see this in action. The $emit_value function performs a set of ref.test checks to see which kind of reference it's dealing with before casting; let's add this line at the very start: It's clearly wrong to assume that $v is a bool reference without first testing it; this is just for demonstration purposes. Without setting any breakpoints, recompiling this code with watgo and reloading the page, we get: The debugger stopped at the instruction causing the exception; moreover, in the Scope pane on the right we can see that the actual type of $v is (ref $Pair) , so it's immediately clear what's going on. I've found this capability extremely valuable when writing (or emitting from a compiler) non-trivial chunks of WASM code using gc types and instructions. "Should I use a debugger or just printfs" is a common topic of debate among programmers. While I'm usually in the "printf debugging" camp, I'm not dogmatic, and will certainly reach for a debugger when the situation calls for it. Specifically, when investigating reference exceptions in WASM, two strong factors tilt the decision towards using a debugger: In general, WASM's printf capabilities aren't great. We can import print-like functions from the host (and - in fact - our sample does just that), but they're not very flexible and dealing with strings in WASM is painful in general. This is compounded even more when working with gc types, because these aren't even visible to the host (they're opaque references). If we want to do printf debugging of gc values, we have to build a lot of scaffolding first. Exception debugging - in general - is much easier with a supportive debugger in hand. Our ref.cast exception from the example above could have happened anywhere in the code. Imagine having to debug a very large WASM program (emitted by a compiler) to find the source of a failed ref.cast ; the debugger takes you right to the spot! In fact, even for C programming, I've always found gdb most useful for pinpointing the source of segmentation faults and similar crashes. In general, WASM's printf capabilities aren't great. We can import print-like functions from the host (and - in fact - our sample does just that), but they're not very flexible and dealing with strings in WASM is painful in general. This is compounded even more when working with gc types, because these aren't even visible to the host (they're opaque references). If we want to do printf debugging of gc values, we have to build a lot of scaffolding first. Exception debugging - in general - is much easier with a supportive debugger in hand. Our ref.cast exception from the example above could have happened anywhere in the code. Imagine having to debug a very large WASM program (emitted by a compiler) to find the source of a failed ref.cast ; the debugger takes you right to the spot! In fact, even for C programming, I've always found gdb most useful for pinpointing the source of segmentation faults and similar crashes.

0 views

watgo - a WebAssembly Toolkit for Go

I'm happy to announce the general availability of watgo - the W eb A ssembly T oolkit for G o. This project is similar to wabt (C++) or wasm-tools (Rust), but in pure, zero-dependency Go. watgo comes with a CLI and a Go API to parse WAT (WebAssembly Text), validate it, and encode it into WASM binaries; it also supports decoding WASM from its binary format. At the center of it all is wasmir - a semantic representation of a WebAssembly module that users can examine (and manipulate). This diagram shows the functionalities provided by watgo: watgo comes with a CLI, which you can install by issuing this command: The CLI aims to be compatible with wasm-tools [1] , and I've already switched my wasm-wat-samples projects to use it; e.g. a command to parse a WAT file, validate it and encode it into binary format: wasmir semantically represents a WASM module with an API that's easy to work with. Here's an example of using watgo to parse a simple WAT program and do some analysis: One important note: the WAT format supports several syntactic niceties that are flattened / canonicalized when lowered to wasmir . For example, all folded instructions are lowered to unfolded ones (linear form), function & type names are resolved to numeric indices, etc. This matches the validation and execution semantics of WASM and its binary representation. These syntactic details are present in watgo in the textformat package (which parses WAT into an AST) and are removed when this is lowered to wasmir . The textformat package is kept internal at this time, but in the future I may consider exposing it publicly - if there's interest. Even though it's still early days for watgo, I'm reasonably confident in its correctness due to a strategy of very heavy testing right from the start. WebAssembly comes with a large official test suite , which is perfect for end-to-end testing of new implementations. The core test suite includes almost 200K lines of WAT files that carry several modules with expected execution semantics and a variety of error scenarios exercised. These live in specially designed .wast files and leverage a custom spec interpreter. watgo hijacks this approach by using the official test suite for its own testing. A custom harness parses .wast files and uses watgo to convert the WAT in them to binary WASM, which is then executed by Node.js [2] ; this harness is a significant effort in itself, but it's very much worth it - the result is excellent testing coverage. watgo passes the entire WASM spec core test suite. Similarly, we leverage wabt's interp test suite which also includes end-to-end tests, using a simpler Node-based harness to test them against watgo. Finally, I maintain a collection of realistic program samples written in WAT in the wasm-wat-samples repository ; these are also used by watgo to test itself. Parse: a parser from WAT to wasmir Validate: uses the official WebAssembly validation semantics to check that the module is well formed and safe Encode: emits wasmir into WASM binary representation Decode: read WASM binary representation into wasmir

0 views

Summary of reading: January - March 2026

"Intellectuals and Society" by Thomas Sowell - a collection of essays in which Sowell criticizes "intellectuals", by which he mostly means left-leaning thinkers and opinions. Interesting, though certainly very biased. This book is from 2009 and focuses mostly on early and mid 20th century; yes, history certainly rhymes. "The Hacker and the State: Cyber Attacks and the New Normal of Geopolitics" by Ben Buchanan - a pretty good overview of some of the the major cyber-attacks done by states in the past 15 years. It doesn't go very deep because it's likely just based on the bits and pieces that leaked to the press; for the same reason, the coverage is probably very partial. Still, it's an interesting and well-researched book overall. "A Primate's Memoir: A Neuroscientist’s Unconventional Life Among the Baboons" by Robert Sapolsky - an account of the author's years spent researching baboons in Kenya. Only about a quarter of the book is really about baboons, though; mostly, it's about the author's adventures in Africa (some of them surely inspired by an intense death wish) and his interaction with the local peoples. I really liked this book overall - it's engaging, educational and funny. Should try more books by this author. "Seeing Like a State" by James C. Scott - the author attempts to link various events in history to discuss "Why do well-intentioned plans for improving the human condition go tragically awry?"; discussing large state plans like scientific forest management, building pre-planned cities and mono-colture agriculture. Some of the chapters are interesting, but overall I'm not sure I'm sold on the thesis. Specifically, the author mixes in private enterprises (like industrial agricultire in the West) with state-driven initiatives in puzzling ways. "Karate-Do: My Way of Life" by Gichin Funakoshi - short autobiography from the founder of modern Shotokan Karate. It's really interesting to find out how recent it all is - prior to WWII, Karate was an obscure art practiced mostly in Okinawa and a bit in other parts of Japan. The author played a critical role in popularizing Karate and spreading it out of Okinawa in the first half of the 20th century. The writing is flowing and succinct - I really liked this book. "A Tale of a Ring" by Ilan Sheinfeld (read in Hebrew) - a multi-generational fictional saga of two families who moved from Danzig (today Gdansk in Poland) to Buenos Aires in late 19th century, with a touch of magic. Didn't like this one very much. "The Wide Wide Sea: Imperial Ambition, First Contact and the Fateful Final Voyage of Captain James Cook" by Hampton Sides - a very interesting account of Captain Cook's last voyage (the one tasked with finding a northwest passage around Canada). The book has a strong focus on his interaction with Polynesian peoples along the way, especially on Hawaii (which he was the first European to visit). "The Suitcase" by Sergei Dovlatov - (read in Russian) a collection of short stories in Dovlatov's typical humorist style. Very nice little book. "The Second Chance Convenience Store" by Kim Ho-Yeon - a collection of connected stories centered around a convenience store in Seoul, and an unusual new employee that began working night shifts there. Short and sweet fiction, I enjoyed it. "A History of the Bible: The Story of the World's Most Influential Book" by John Barton - a very detailed history of the Bible, covering both the old and new testaments in many aspects. Some parts of the book are quite tedious; it's not an easy read. Even though the author tries to maintain a very objective and scientific approach, it's apparent (at least for an atheist) that he skirts as close as possible to declaring it all nonsense, given that he's a priest! "Rust Atomics and Locks: Low-Level Concurrency in Practice" by Mara Bos - an overview of low-level concurrency topics using Rust. It's a decent book for people not too familiar with the subject; I personally didn't find it too captivating, but I do see the possibility of referring to it in the future if I get to do some lower-level Rust hacking. A comment on the code samples: it would be nice if the accompanying repository had test harnesses to observe how the code behaves, and some benchmarks. Without this, many claims made in the book feel empty without real data to back them up, and it's challenging to play with the code and see it perform in real life. "Hot Chocolate on Thursday" by Michiko Aoyama - a bit similar to "What You Are Looking for Is in the Library" by the same author: connected short stories about ordinary people living their life in Japan (with one detour to Australia). Slightly worse than the previous book, but still pretty good. "The Silmarillion" by J.R.R. Tolkien - enen though I'm a big LOTR fan, I've never gotten myself to read this one, due to its reputation for being difficult. What changed things eventually (25 years after my first read through of LOTR) is my kids! They liked LOTR so much that they went straight ahead to Silmarillion and burned through it as well, so I couldn't stay behind. What can I say, this book is pretty amazing. The amazing thing is how a book can be both epic and borderline unreadable at the same time :) Tolkien really let himself go with the names here (3-4 new names introduced per page, on average), names for characters, names for natural features like forests and rivers, names for all kinds of magical paraphenalia; names that change in time, different names given to the same thing by different peoples, and on and on. The edition I was reading has a helpful name index at the end (42 pages long!) which was very helpful, but it still made the task only marginally easier. Names aside though, the book is undoubtedly monumental; the language is outstanding. It's a whole new mythology, Bible-like in scope, all somehow more-or-less consistent (if you remember who is who, of course); it's an injustice to see this just as a prelude to the LOTR books. Compared to the scope of the Simlarillion, LOTR is just a small speck of a quest told in detail; The Silmarillion - among other things - includes brief tellings of at least a dozen stories of similar scope. Many modern book (or TV) series build whole "universes" with their own rules, history and aesthetic. The Silmarillion must be considered the OG of this. "Travels with Charley in Search of America" by John Steinbeck "Deep Work" by Cal Newport "The Philadelphia chromosome" by Jessica Wapner "The Price of Privelege" by Madeline Levine

0 views

Notes on Lagrange Interpolating Polynomials

Polynomial interpolation is a method of finding a polynomial function that fits a given set of data perfectly. More concretely, suppose we have a set of n+1 distinct points [1] : And we want to find the polynomial coefficients {a_0\cdots a_n} such that: Fits all our points; that is p(x_0)=y_0 , p(x_1)=y_1 etc. This post discusses a common approach to solving this problem, and also shows why such a polynomial exists and is unique. When we assign all points (x_i, y_i) into the generic polynomial p(x) , we get: We want to solve for the coefficients a_i . This is a linear system of equations that can be represented by the following matrix equation: The matrix on the left is called the Vandermonde matrix . This matrix is known to be invertible (see Appendix for a proof); therefore, this system of equations has a single solution that can be calculated by inverting the matrix. In practice, however, the Vandermonde matrix is often numerically ill-conditioned, so inverting it isn’t the best way to calculate exact polynomial coefficients. Several better methods exist. Lagrange interpolation polynomials emerge from a simple, yet powerful idea. Let’s define the Lagrange basis functions l_i(x) ( i \in [0, n] ) as follows, given our points (x_i, y_i) : In words, l_i(x) is constrained to 1 at and to 0 at all other x_j . We don’t care about its value at any other point. The linear combination: is then a valid interpolating polynomial for our set of n+1 points, because it’s equal to at each (take a moment to convince yourself this is true). How do we find l_i(x) ? The key insight comes from studying the following function: This function has terms (x-x_j) for all j\neq i . It should be easy to see that l'_i(x) is 0 at all x_j when j\neq i . What about its value at , though? We can just assign into l'_i(x) to get: And then normalize l'_i(x) , dividing it by this (constant) value. We get the Lagrange basis function l_i(x) : Let’s use a concrete example to visualize this. Suppose we have the following set of points we want to interpolate: (1,4), (2,2), (3,3) . We can calculate l'_0(x) , l'_1(x) and l'_2(x) , and get the following: Note where each l'_i(x) intersects the axis. These functions have the right values at all x_{j\neq i} . If we normalize them to obtain l_i(x) , we get these functions: Note that each polynomial is 1 at the appropriate and 0 at all the other x_{j\neq i} , as required. With these l_i(x) , we can now plot the interpolating polynomial p(x)=\sum_{i=0}^{n}y_i l_i(x) , which fits our set of input points: We’ve just seen that the linear combination of Lagrange basis functions: is a valid interpolating polynomial for a set of n+1 distinct points (x_i, y_i) . What is its degree? Since the degree of each l_i(x) is , then the degree of p(x) is at most . We’ve just derived the first part of the Polynomial interpolation theorem : Polynomial interpolation theorem : for any n+1 data points (x_0,y_0), (x_1, y_1)\cdots(x_n, y_n) \in \mathbb{R}^2 where no two x_j are the same, there exists a unique polynomial p(x) of degree at most that interpolates these points. We’ve demonstrated existence and degree, but not yet uniqueness . So let’s turn to that. We know that p(x) interpolates all n+1 points, and its degree is . Suppose there’s another such polynomial q(x) . Let’s construct: That do we know about r(x) ? First of all, its value is 0 at all our , so it has n+1 roots . Second, we also know that its degree is at most (because it’s the difference of two polynomials of such degree). These two facts are a contradiction. No non-zero polynomial of degree \leq n can have n+1 roots (a basic algebraic fact related to the Fundamental theorem of algebra ). So r(x) must be the zero polynomial; in other words, our p(x) is unique \blacksquare . Note the implication of uniqueness here: given our set of n+1 distinct points, there’s only one polynomial of degree \leq n that interpolates it. We can find its coefficients by inverting the Vandermonde matrix, by using Lagrange basis functions, or any other method [2] . The set P_n(\mathbb{R}) consists of all real polynomials of degree \leq n . This set - along with addition of polynomials and scalar multiplication - forms a vector space . We called l_i(x) the "Lagrange basis" previously, and they do - in fact - form an actual linear algebra basis for this vector space. To prove this claim, we need to show that Lagrange polynomials are linearly independent and that they span the space. Linear independence : we have to show that implies a_i=0 \quad \forall i . Recall that l_i(x) is 1 at , while all other l_j(x) are 0 at that point. Therefore, evaluating s(x) at , we get: Similarly, we can show that a_i is 0, for all \blacksquare . Span : we’ve already demonstrated that the linear combination of l_i(x) : is a valid interpolating polynomial for any set of n+1 distinct points. Using the polynomial interpolation theorem , this is the unique polynomial interpolating this set of points. In other words, for every q(x)\in P_n(\mathbb{R}) , we can identify any set of n+1 distinct points it passes through, and then use the technique described in this post to find the coefficients of q(x) in the Lagrange basis. Therefore, the set l_i(x) spans the vector space \blacksquare . Previously we’ve seen how to use the \{1, x, x^2, \dots x^n\} basis to write down a system of linear equations that helps us find the interpolating polynomial. This results in the Vandermonde matrix . Using the Lagrange basis, we can get a much nicer matrix representation of the interpolation equations. Recall that our general polynomial using the Lagrange basis is: Let’s build a system of equations for each of the n+1 points (x_i,y_i) . For : By definition of the Lagrange basis functions, all l_i(x_0) where i\neq 0 are 0, while l_0(x_0) is 1. So this simplifies to: But the value at node is , so we’ve just found that a_0=y_0 . We can produce similar equations for the other nodes as well, p(x_1)=a_1 , etc. In matrix form: We get the identity matrix; this is another way to trivially show that a_0=y_0 , a_1=y_1 and so on. Given some numbers \{x_0 \dots x_n\} a matrix of this form: Is called the Vandermonde matrix. What’s special about a Vandermonde matrix is that we know it’s invertible when are distinct. This is because its determinant is known to be non-zero . Moreover, its determinant is [3] : Here’s why. To get some intuition, let’s consider some small-rank Vandermonde matrices. Starting with a 2-by-2: Let’s try 3-by-3 now: We can use the standard way of calculating determinants to expand from the first row: Using some algebraic manipulation, it’s easy to show this is equivalent to: For the full proof, let’s look at the generalized n+1 -by- n+1 matrix again: Recall that subtracting a multiple of one column from another doesn’t change a matrix’s determinant. For each column k>1 , we’ll subtract the value of column k-1 multiplied by from it (this is done on all columns simultaneously). The idea is to make the first row all zeros after the very first element: Now we factor out x_1-x_0 from the second row (after the first element), x_2-x_0 from the third row and so on, to get: Imagine we erase the first row and first column of . We’ll call the resulting matrix . Because the first row of is all zeros except the first element, we have: Note that the first row of has a common factor of x_1-x_0 , so when calculating \det(W) , we can move this common factor out. Same for the common factor x_2-x_0 of the second row, and so on. Overall, we can write: But the smaller matrix is just the Vandermonde matrix for \{x_0 \dots x_{n-1}\} . If we continue this process by induction, we’ll get: If you’re interested, the Wikipedia page for the Vandermonde matrix has a couple of additional proofs.

0 views