Results for 2008

Computational physics is one branch of physics that solves problems in physics by implementing computational algorithms. Especially in theoretical physics, the study of numerical physics is significantly useful and often used for simulations.

One of the most exiting branches of computational physics is numerical relativity. It is the astrophysical study to simulate objects in space-time universe like stars, black holes, and gravitational waves, based on Einstein’s Theory of General Relativity.

The physical law in space-time can basically be represented by one equation: the Einstein Equation. Therefore, numerical relativity is primarily the study of numerical solution of the Einstein equations:

Gij = 8πTij.

These are16 coupled hyperbolic-elliptic nonlinear partial differential equations (which reduce to 10 by symmetry). They equate the curvature of space-time (Einstein tensor G) with the energy and momentum in the space-time (stress-energy tensor T). As the result, the differential equations evolve time and space derivatives. Simulating the movements of objects in space-time on computer requires the methods to take spatial derivatives and time derivative. I introduce some of the practical methods we learned in class to deal with these problems.

Special derivatives
Einstein equations by themselves cannot be differentiated exactly. Therefore, we need to change them to the forms so that they can be differentiated exactly. Fourier series expansion, expansion in terms of Chebyshev polynomials are used for this purpose. These methods are effective because they have high convergence rates.

Time derivatives
Runge-Kutta equations of fourth-order are often used in this step. As we increase the order of Runge-Kutta equations, we can acquire arbitrary convergence rate.

Numerical Relativity group at the Albert Einstein Institute is expert in these simulations, and I site their resulted images.

General information of Numerical Relativity

Image and Movie Archive of the AEI Numerical Relativity Group

In our class we have talked about single and double precision mathematics floating point numbers and how these are stored. Another class of computer numbers also exist called arbitrary precision numbers. These are unique in that they allow the user to specify how many digits of precision they would like stored.

Regular single and doubles precision numbers are allocated a fixed amount of space in memory, 32 and 64 bits respectively and hence are limited to about 8 or 16 decimal digits of precision. Unlimited or infinite precision numbers on the other hand occupy different amounts of memory depending on the number being stored. These terms however are misnomers since the number of digits that can be stored, while very large is always finite and limited by software and memory constraints.

Arbitrary precision numbers are useful in many specialized tasks from scientific computing to cryptography. The cryptography used for securing web transactions for example often requires working with thousand plus digit numbers. In mathematics they are used when calculating the value of constants such at pi or the square root of 3. Sometimes, to see the impact of rounding errors, it may be useful to perform the same calculation in standard precision and again in extended precision. Drawing fractals in high magnification is another task that benefit from higher precision numbers. Applications where the output is viewed by humans and speed is not a concern, extended precision may also make sense.

Extended precision numbers can mitigate the many types of failures that can result from bad numerical programming. Disasters written about extensively in earlier posts such as failing rockets may have been averted by the usage of higher precision numbers. Big integers for example can prevent overflow errors. Big floating point numbers mitigate the impact of catastrophic cancellation that results from subtracting two similar numbers and can sometimes help address the issue of range reduction.

Bignum classes are built into many progra    mming languages. Java for example has the BigInt and BigFloat classes. specific classes exist in many languges to For example when calculating pi to many decimal places, regular single and double precision numbers are not sufficient.

There is however a cost to the use of arbitrary precision numbers and hence their limited usage. Computer hardware is optimized to store and operate on either single precision or double precision numbers. For example, graphics cards are designed to perform operations on arrays of single precision numbers since pixel colors can be fully represented with only 32 bits.

Thus working with extended precision numbers necessitates the use of software to perform simple operations which can be much slower. Efficient algorithms have been developed to perform common mathematical operations. Addition and subtraction are relatively easy and can be accomplished in O(n) time where n is the number of digits of precision involved. Multiplication is surprising complicated and many algorithms with different time complexity exist. Algorithms with complexity from O(n²) to O(n log n 2^{log* n}) exist.

http://en.wikipedia.org/wiki/Bignum

http://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations

With so many hot applications like computer graphics, advanced audio synthesis, and video relying heavily on matrix multiplication, Apple is making efforts to make this process more efficient and accurate with a dedicated vector processing unit (VPU) that can be carried out these operations for various applications.  Rather than having to simplify the complexity of these matrices to a close approximation, this technology allows for calculations to carried out on the original matrices themselves, without worrying about rounding and truncation errors.

In Apple’s background description of the technology, developed by Ali Sazegari, a description of the main techniques is described:

“The invention provides means for storing in first vector registers within a vector processing system, multiple copies of each value of a first matrix, and storing in second vector registers within a vector processing system, values from each row of a second matrix. Upon storing the values in these first and second vector registers, the present invention dot multiplies each of the first vector registers having multiple copies of each value associated with a first row of the first matrix by each of the second vector registers. The dot product values obtained by this dot multiplication are each added to values in vector registers comprising a product matrix, and the multiplication and addition is again repeated for the first vector registers having values from every row of the first matrix. This method is bit-by-bit compatible with the traditional matrix multiplication, and does not introduce any discrepancies in rounding or otherwise which often occur because of changes in the order of operations.”

This process is carried out by divided very complex matrices, into smaller and easy to compute submatrices. Combined with the VPUs impressive ability to carry  out operations on an entire vector register simultaneously.  If this technology proves successful, it may allow for developers to no longer shy away from using complex matrices for highly accurate information processing.  Despite being serving a very specific and seemingly narrow purpose, matrix multiplication is one of the most computational demanding tasks widely used. If developers didn’t have to put thought into handling these matrices efficiently, then more system resources can be devoted to enhancing other processes or by simply saving performance since the CPU doesn’t have to worry about the computations.

 Sources:

http://www.macsimumnews.com/index.php/archive/apple_patent_is_for_matrix_multiplication_in_a_vector_processing_unit

The National Science Teachers’ Association’s annual national conference, scheduled this week in Boston, MA, will draw thousands of science educators from across the country. Several Beyond Penguins and Polar Bears project staff members are attending the conference. We’ll post any polar-related news or sessions of interest. Check back often if you are unable to attend the conference, or just missed a particular session!

Additionally, we are scheduled to present an IPY Session: Arctic and Antarctic Biology on Sunday, March 30 from 8:00-9:00 a.m. in the Boston Convention and Exhibition Center, Room 257 B. In our presentation, we’ll provide an overview of current news relating to the living systems of the polar regions and a list of news providers that offer RSS feeds and email updates. We’ll also discuss the learning progression found in the National Science Education Standards Life Science Content Standard. We’ll highlight the interactive Science Literacy maps from the National Science Digital Library and provide a resource list of lesson plans, activities, and web sites for K-12 educators.

If you aren’t able to attend, you can download our presentation and resource list here!

Finally, we’ll be at the National Science Digital Library booth (#1145) in the exhibit hall throughout the conference. Stop by to learn more about Beyond Penguins or any of the wonderful resources available from NSDL!

Stay tuned for more news from Boston!

There are many ways to compute the n-th Fibonacci number easily: there is a well-known close formula. Another approach is matrix multiplication.

Let’s take a look at this approach first.

As we all know, F_k+2 = F_k+1 + F_k. Thus, we can expand this equation to a matrix equation:

[1]

Now, it’s easy to see the following:

[2]

Now, we only need to compute this n-th power of a magic matrix for Fibonacci numbers.

If one 2-by-2-matrix-multiplication costs you $1, you only need to pay $n-1 to compute the n-th power. Or, even better, you only need about $log₂ n with some fancy optimization; I believe I learned this in CS 100/211. See the following pseudo-code:

int power(int x, int n) // computes x ^ n
{
if (n == 0) return 1;
if (n % 2 == 1) return power(x, n-1) * x;
int tmp = power(x, n/2);
return tmp * tmp;
}

Then, how about the SUM of first n terms of Fibonacci numbers? How about the sum of numbers between n-th term and m-th term? One easy way is to compute each number and add them up; this will cost you too much.

Let us call the magic matrix above M.

In order to get the sum of first n Fibonacci numbers, we need to compute M, M², M³, … , Mⁿ. This will cost you roughly $n²/2. Can we do better?

We want our final matrix to be M + M² + M³ + … + Mⁿ. Now recall how we compute xⁿ with O(log n) multiplications. Suppose we have the matrix M + M² + M³ + … + M^n/2. If we multiply M^n/2 to this sum of matrices, and then add it to the original sum, we finally get M + M² + M³ + … + Mⁿ. We are using similar optimization here.

Matrix powerSum(M, n)
{
if (n == 0) return E; // unit 2-by-2 matrix
if (n == 1) return M;
if (n % 2 == 1)
{
S = powerSum(M, (n-1)/2);
T = S * M^(n-1)/2; // we can optimize computing M^(n-1)/2, too.
return S + T + M^n; // M^n can be computed from M^(n-1)/2.

} else {
S = powerSum(M, n/2);
T = S * M^(n/2);
return S + T;
}
}

You can easily see that we do NOT use matrix multiplication n^2 times; it’s O(log n) if we optimize computation of M^(n/2) - we can do this by forcing this powerSum method to return M^n as well when it returns the sum of M, M^2, …, M^n - then, the caller can use this result to compute M^2n or M^(2n+1), and so on.

References

[1] http://upload.wikimedia.org/math/3/d/c/3dcf66f1eb82acf058a19800b2094fa9.png

[2] http://upload.wikimedia.org/math/a/6/0/a6083f85f39b468210f5715a8e30d572.png

What can we do with ODEs? It turns out the ODEs are an integral part of the cloth simulation process.

For simulation purposes cloth is often modeled as a quad mesh whose vertices are connected by stiff springs. In general cloth resists stretching, this property is called extensibility. In order to simulate this, the springs must be made stiffer and stiffer in the model, and as stiffness increases numerical methods for simulating the springs in cloth become intractable [1].

Another approach builds on this by adding constrained dynamics to enforce inextensibility of cloth. Constraints include edge length preservation in the cloth. This approach incorporates two topics from scientific computing, ordinary differential equation solvers and root finding methods. The general process can be described as follows [2]:

1) Start with some initial position of the cloth x0 and the world model.
2) Solve all constraints simultaneously.
3) Apply collision response.
4) Implicitly integrate the springs.

While there is some somewhat fancier math involved, solving the implicit constraint dynamics in step 2 boils down to finding the zeros of an equations map calculated points to points which preserve the length of edges of the cloth. As we know from class, root finding can be performed by numerous methods and in this case Newton’s is generally used.
Additionally, the cloth simulation is performed in discrete time steps. In this setup in order to perform the implicit integration of the springs and ODE solver is needed. An equation of the form below must be solved.

x’’ = (M-1) * f(x, x’)

As is popular in many areas of engineering and simulation, RK4 is the method of choice due to its speed to stability ratio.

The result of each of these steps is life-like cloth simulation with realistic folds and creases that maintains its length without stretching unnecessarily.

[1] Kwang-Jin Choi and Hyeong-Seok Ko. Stable but responsive cloth. In John Hughes, editor, SIGGRAPH 2002 Conference Proceedings, Annual Conference Series, pages 604–611. ACM Press/ACM SIGGRAPH, 2002.
[2] Goldenthal, R., Harmon, D., Fattal, R., Bercovier, M., and Grinspun, E. 2007. Efficient simulation of inextensible cloth. ACM Trans. Graph. 26, 3 (Jul. 2007), 49.

The singular value decomposition we learned in class can be used as a way to compress images. If you were given a 256 by 256 image, you could represent it by storing RGB values in a 256 by 256 matrix which correspond to each pixel in the image. In order to transmit this image without any compression, you would need to transmit 65536 RGB values.

In order to compress this image, it is possible to use singular value decomposition. The first thing to do is calculate the SVD as discussed in class of the matrix that represents the image. This will give something of the form:

A=UΣV*

Where A is the original m by n matrix (in the case above 256 by 256), U is an m by m matrix containing an orthonormal basis of an m dimensional space, Σ is a diagonal matrix containing the singular values of A in decreasing order, and V* is an n by n matrix with columns that form an orthonormal basis of an n dimensional space.

Because the singular values of most matrices become very small after the first few, it is possible to simply not include them in the SVD in order to make there be less information to transmit. When this is done, the corresponding parts of U and V* are also not included. If only 50 singular values were used for the image above, only 50 vectors from U and 50 vectors from V* would need to be transferred. Since each vector in U and V* has 256 values in them, this can get rid of a lot of values, while not losing much information of the image. A more detailed explanation of how this works with some images showing different amounts of singular values used is shown at the University of Wisconsin La Crosse website below. There is also a link to a paper from the University of Calgary about using SVD for image compression that explores some more advanced methods than I described.

http://www.uwlax.edu/faculty/will/svd/compression/index.html
http://math.ucalgary.ca/~laf/teaching/Material/Arnold.pdf

“NUMBERS of the form nsqrt{-1} are imaginary, but can still be used in equations.”^[1]
After seeing this, I was really stunned at what in the world that equation means.
Later, I found another weird equations through xkcd #217^[2].
I was really confused and I decided to find out what these equations mean and/or how they can be computed, really.
Especially, “e to the pi times i” was really hard to interpret.
I think a few of you already know about this thing, but for those like me, this might be really interesting.
The way ‘people’ calculate the value for “e to the pi times i” is… involving Taylor’s Series, which we are using these days in class.
In short, by Taylor’s Series, we can do the following by plugging in xi in place of x:

^[3]

The very last right side of the equation tells us that e^{pi i} = cospi + i sinpi = -1.
So, whatever it means that e has an imaginary power, the result is… an integer number!

If you find out this very interesting, there is another explanation to this using Calculus - see reference[3].

[1] http://xkcd.com/179/

[2] http://xkcd.com/217/

[3] http://www.math.toronto.edu/mathnet/questionCorner/epii.html

For our first project we have used ray tracing to create images of objects. As a result we were able to get decent images of long extinct water(?) bunnies. However, aside from the root calculation, tracing an object took a long time. It was obvious that once the images became complex, ray tracing would take a much longer time. For many years ray tracing was limited to offline-rendering, since the computations took so much time. Therefore real-time 3D computer graphics had come up with an alternate solution: Rasterization

As the most popular method of rendering 3D images on a computer, rasterization is a method of taking an image described in vector graphics format and converting it to raster graphic images (e.g. SVG to BMP). This meant just converting the 3D geometry into 2D pixels, and lighting/shading had to be taken care of other methods. Because it was simple, rasterization was very fast compared to ray tracing, and thus well suited to the demands of various applications. The development of graphic cards and DirectX, OpenGL allowed rasterization to become the most common method used in image rendering.

However, with the development of OpenRT Ray tracing library (Saarland University), real-time ray tracing was possible. (Well documented in other posts) The author Daniel Pohl re-introduced ray tracing into the well known 3D game Quake 3 & 4 to see its superiority over rasterization. Through his works of various pictures and videos, we can see how rasterzation can create lots of graphic/lighting artifacts, especially noticeable in shadows, while ray tracing technology can achieve outstanding, realistic images through a few simple lines of code.

The possibly biggest advantage that ray tracing brings in the era of multicore/cluster CPU is that it is perfect for parallelization.
Calculations for each ray are independent, and can be done in parallel. As processors are moving toward many small CPUs, ray tracing will be perfect for future computers. His works show that ray-tracing has parallel efficiency of almost 1.

Companies like id Software put in a lot of effort in optimizing rendering engines so that the images look realistic as possible.
However, even with the latest technologies such as megatexture and the id5 tech engine, it is clear in the pictures that ray tracing is superior in rendering 3D images in many ways. Although there are problems yet to be overcome, ray tracing will sooner or later take over as the main method for 3D graphic rendering. It is possible that ray tracing will lead GPUs into multicore to maximize parallelization.

Sources:
Ray Tracing and Gaming - Quake 4: Ray Traced Project
Ray Tracing and Gaming - One Year Later

Taking a hint from rock stars, the researchers from Polar-Palooza, a multimedia project supported by the National Science Foundation and the National Aeronautics and Space Administration, are going on national tour of some 14 cities this year.

In three-day appearances at each location, the “cast” of scientists will present “Stories from a Changing Planet,” with video and authentic props, workshops for educators, and briefings for the press. The tour starts off in Washington, D.C.,  on March 13-14 at the National Geographic Museum and crisscrosses the country into the fall. Find the 2008 Venues at http://passporttoknowledge.com/polar-palooza/pp04.php.

Professor Philipp Slusallek at the University of Saarland, and other members of the OpenRT project have developed algorithms that have drastically sped up how fast ray tracing can be done. Using these algorithms on several computers (equivalent to about 36 GHz) the OpenRT group was able to achieve 20 fps while running a ray traced version of Quake 3. While it would still be impossible to run something like this on any normal computer, it is a considerable step towards being able to use ray tracing in gaming and other computer graphic applications.

The OpenRT project does not only focus on making more efficient ray tracing algorithms. They also have done a lot of work on making the finished image be as realistic as possible. Their algorithms even account for complicated light behavior, such as caustics and multiple reflections. Images that demonstrate what OpenRT can do can be viewed on their website, which is listed below.

OpenRT is relevant to this class because our first project is about ray tracing. Also, ray tracing is essentially finding the roots of functions, which is what we have focused on so far.

OpenRT’s homepage: http://www.openrt.de/
BBC story on OpenRT:http://news.bbc.co.uk/2/hi/technology/6457951.stm

Helen, a very clever girl, never makes any mistake on her homework. Her GPA used to be 4.3 until she got one miserable A due to a silly numerical computing error… it was a disaster.. a nightmare.. even now, she sometimes gets pissed off when she thinks about it. So, what happened to Helen last semester?

Last semester, she was working on her final project report on MS Excel 2007. The project report contained a lot of numerical values, both integers and rationals, so that she couldn’t check all the values manually. The first data set she got from her experiment was following:

So far, it seems fine (really?). The data contained more than 200 rows so that she couldn’t validate if the values are correct or not. Later on, using this data, she calculated something more complicated and then finally made a super-awesome graph! Helen was so satisfied with her work because she spent 5 days to get the most accurate data she could achieve and then spent another whole day to write up this beautiful report with this awesome graph!

On the following day, she was presenting her final report in class. Well, in the middle of her presentation, somebody raised his hand up and asked,

“Hey, how did you get that number? You got 77.1 * 850 = 100,000, but I don’t think it’s correct. It should be less than 77,000 because 77.1 * 1,000 is only 77,000!”

Confused and puzzled, Helen said,

“Oh, well, I don’t know… Um… I just used MS Excel to do the calculations… I didn’t do it so I think it’s MS Excel’s bug..?”^{[1], [2]}

Anyway, her presentation was not good at all and she did not get a good grade on her final project: ‘due to a careless mistake when validating the data’. Even if she did so well on the exams and other projects, she did not receive A+ because of her final presentation and errors on the data. Now her GPA is 4.28 and she is so sad.

Well, why is this happening? That’s because YOUR COMPUTER SOMETIMES DOES NOT REALLY KNOW HOW TO EXPRESS FLOATING NUMBERS IN BINARY NUMBERS! If you think about all the famous disasters caused by numerical errors, many of them can be categorized in one group: errors when dealing with floating numbers. If you check out the websites [1] and [2], you will find out why this bug was made (even on Microsoft Excel 2007, 2007! MS!).

Excel:

1.to surpass others or be superior in some respect or area; do extremely well: to excel in math.^[3]

Personally, I hate to use float/double (in C++) kinds of primitive types because of this imperfection. Whenever I calculate something like 1/3 * 3, the result is not 1, but 0.999999999999999! This sometimes really bugs me down so that I have a weird habit now: try to avoid using any float-type data (what? how?). For instance, if you calculate Euclidean distance between two points, you will write something like this:

double dist = sqrt ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));

If x1, x2, y1, and y2 are integers, I would rather write

int sqr_dist = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);

and keep this squared distance until I really need to calculate the actual value.

Another thing is, I always use some small constant ‘epsilon’ when I compare two non-integer values:

const double eps = 1.e-6; // 10^-6

bool same(double A, double B) { return abs(A-B) < eps; }

I never write something like:

if (A == B) // …

if A and B are not integer-type because 1.0 is different from 0.999999999999…!

After all, my point is, whoever realized the existence of irrational numbers^[4], they made our computers and engineers suffer.. (since we can represent any rational as a pair of two integers, rationals are OK).

P.S. I hope this first posting can be a good reason for the professor to accept my late homework #1…?

[Sources]

[1] http://www.joelonsoftware.com/items/2007/09/26b.html

[2] http://blogs.msdn.com/excel/archive/2007/09/25/calculation-issue-update.aspx

[3] http://dictionary.reference.com/browse/excel

[4] http://en.wikipedia.org/wiki/Real_number#History