Month: February 2012

Dot, cross, and quaternion products

Posted on by John

This post will show that

quaternion product = cross product – dot product.

First, I’ll explain what quaternions are, then I’ll explain what the equation above means.

The complex numbers are formed by adding to the real numbers a special symbol i with the rule that i2 = -1. The quaternions are similarly formed by adding to the real numbers i, j, and k with the requirement [1] that

i2 = j2 = k2 = ijk = -1.

A quaternion is a formal sum of a real number and real multiples of the symbols i, j, and k. For example,

q = w + xi + yj + zk.

In this example we say w is the real part of q and xi + yj + zk is the vector part of q, analogous to the real and imaginary parts of a complex number.

The quaternion product of two vectors (x, y, z) and (x´, y ´, z´) is the product of q = xi + yj + zk and q‘ = x’i + y’j + z’k as quaternions. The quaternion product qq´ works out to be

– (xx´ + yy´ + zz´) + (yz´ – zy´)i +(zx´ – xz´)j + (xy´ – yx´)k

The real part is the negative of the dot product of (x, y, z) and (x´, y´, z´) and the vector part is the cross product of (x, y, z) and (x´, y´, z´).

NB: The context here is vectors of length 3 embedded in the quaternions. In general the quaternion product is more complicated. The products are simpler in the special case of the real component being zero. See, for example, this post.

This relationship is an interesting bit of algebra on its own, but it is also historically important. In the 19th century, a debate raged regarding whether quaternions or vectors were the best way to represent things such as electric and magnetic fields. The identity given here shows how the two approaches are related.

Related posts

[1] To multiply two quaternions, you need to know how to multiply i, j, and k by each other. It is not immediately obvious, but you can derive everything from i2 = j2 = k2 =ijk = -1. For example, start with ijk = -1 and multiply both sides on the right by k. So ijk2 = –k, and since k2 = -1, ij = k. Similar manipulations show jk = i and ki = j.

Next, (jk)(ki) = ij, but it also equals jk2i = –ji, so ij = –ji = k. Similarly kj = –jk = –i and ik = –ki = –j and this completes the multiplication table for i, j, and k.

The way to remember these products is to imagine a cycle: i -> j -> k -> i. The product two consecutive symbols is the next symbol in the cycle. And when you traverse the cycle backward, you change the sign.

Irreproducible research on 60 Minutes

Posted on by John

If your research cannot be reproduced, you might end up on 60 Minutes. Two days ago the new show ran a story about irreproducible research at Duke. You can find the video clip here.

I believe the 60 Minutes piece was somewhat misleading. It focused on data manipulation and implied that the controversial results followed from the manipulated data. As Keith Baggerly explains here, that is not the case. The conclusions do not follow from the (erroneous) data. The analysis itself was irreproducible. That discovery started the whole saga.

Update: Here’s some footage that 60 Minutes recorded but did not include on Sunday. “The systems we have in academia, especially with something this complicated, shield sloppy science and fraud.”

Update: Guess someone took that video down. Sorry.

More posts on reproducibility

Hard science, soft science, hardware, software

Posted on by John

The hard sciences—physics, chemistry, astronomy, etc.—boasted remarkable achievements in the 20th century. The credibility and prestige of all science went up as a result. Academic disciplines outside the sciences rushed to append “science” to their names to share in the glory.

Science has an image of infallibility based on the success of the hard sciences. When someone says “You can’t argue with science,” I’d rather they said “It’s difficult to argue with hard science.”

The soft sciences get things wrong more often. Sciences such as biology and epidemiology — soft compared to physics, but hard compared to sociology — often get things wrong. In softer sciences, research results might be not even wrong.

I’m not saying that the softer sciences are not valuable; they certainly are. Nor am I saying they’re easier; in some sense they’re harder than the so-called hard sciences. The soft sciences are hard in the sense of being difficult, but not hard in the sense of studying indisputably measurable effects and making sharp quantitative predictions. I am saying that the soft sciences do not deserve the presumption of certainty they enjoy by association with the hard sciences.

There’s a similar phenomena in computing. Computing hardware has made astonishing progress. Software has not, but it enjoys some perception of progress by association. Software development has improved over the last 60 years, but has made nowhere near the progress of hardware (with a few exceptions). Software development has gotten easier more than it has gotten better. (Old tasks have gotten easier to do, but software is expected to do new things, so it’s debatable whether all told software development has gotten easier or harder.)

Related posts

Technology is more than computers

Posted on by John

Moira Gunn recently interviewed Michel Stipe and Mike Mills of R.E.M. about how they have incorporated technology into their music.

Well there really wasn’t a lot of technology to be thought about when we started. You plugged your amps in and you plugged your guitars into that and hopefully your microphone worked and that was it.

R.E.M. was making music just like our stone age ancestors with only amplifiers, electric guitars, and microphones.

Alan Kay defined technology as “anything that was invented after you were born.” Danny Hillis defined it as “everything that doesn’t work yet.” So if you were born after the invention of amplifiers, electric guitars, and microphones, these reliable technologies aren’t technology.

Why the symbol for magnetic field is ‘B’

Posted on by John

I asked on Twitter the other day

What is the historical reason for denoting magnetic filed “B”?

Eric Eekhoff sent me an answer and with his permission I’m copying his email below:

Hi John,

I saw your question on your GrokEM twitter account about why magnetic field was denoted as B. I recall that Maxwell just used the letters A through H for vectors in his Treatise on Electricity and Magnetism and some of them stuck and some of them didn’t. A is still used for vector potential, B for magnetic field (or magnetic induction or flux density, depending who you ask), H for magnetic intensity, etc. Maxwell used C and G for other vectors that I don’t recall at the moment. They, for some reason, never stuck.

Hope that helps. Have a good day,

Eric

Mixing R, Python, and Perl in 14 lines of code

Posted on by John

This is a continuation of my previous post, Running Python and R inside Emacs. That post shows how to execute independent code blocks in Emacs org-mode. This post illustrates calling one code block from another, each written in a different language.

The example below computes sin2(x) + cos2(x) by computing the sine function in R, the cosine function in Python, and summing their squares in Perl. As you’d hope, it returns 1. (Actually, it returns 0.99999999999985 on my machine.)

To execute the code, go to the #+call line and type C-c C-c.

#+name: sin_r(x=1)
#+begin_src R
sin(x)
#+end_src

#+name: cos_p(x=0)
#+begin_src python
import math
return math.cos(x)
#+end_src

#+name: sum_sq(a = 0, b = 0)
#+begin_src perl
$a*$a + $b*$b;
#+end_src

#+call: sum_sq(sin_r(1), cos_p(1))

Apparently each function argument has to have a default value. If that’s documented, I missed it. I gave the sine and cosine functions default values that would cause the call to sum_sq to return more than 1 if the defaults were used.

Running Python and R inside Emacs

Posted on by John

Emacs org-mode lets you manage blocks of source code inside a text file. You can execute these blocks and have the output display in your text file. Or you could export the file, say to HTML or PDF, and show the code and/or the results of executing the code.

Here I’ll show some of the most basic possibilities. For much more information, see  orgmode.org. And for the use of org-mode in research, see A Multi-Language Computing Environment for Literate Programming and Reproducible Research.

Source code blocks go between lines of the form

    #+begin_src
    #+end_src

On the #+begin_src line, specify the programming language. Here I’ll demonstrate Python and R, but org-mode currently supports C++, Java, Perl, etc. for a total of 35 languages.

Suppose we want to compute √42 using R.

    #+begin_src R
    sqrt(42)
    #+end_src

If we put the cursor somewhere in the code block and type C-c C-c, org-mode will add these lines:

    #+results:
    : 6.48074069840786

Now suppose we do the same with Python:

    #+begin_src python
    from math import sqrt
    sqrt(42)
    #+end_src

This time we get disappointing results:

    #+results:
    : None

What happened? The org-mode manual explains:

… code should be written as if it were the body of such a function. In particular, note that Python does not automatically return a value from a function unless a return statement is present, and so a ‘return’ statement will usually be required in Python.

If we change sqrt(42) to return sqrt(42) then we get the same result that we got when using R.

By default, evaluating a block of code returns a single result. If you want to see the output as if you were interactively using Python from the REPL, you can add :results output :session following the language name.

    #+begin_src python :results output :session
    print "There are %d hours in a week." % (7*24)
    2**10
    #+end_src

This produces the lines

    #+results:
    : There are 168 hours in a week.
    : 1024

Without the :session tag, the second line would not appear because there was no print statement.

I had to do a couple things before I could get the examples above to work. First, I had to upgrade org-mode. The version of org-mode that shipped with Emacs 23.3 was quite out of date. Second, the only language you can run by default is Emacs Lisp. You have to turn on support for other languages in your .emacs file. Here’s the code to turn on support for Python and R.

    (org-babel-do-load-languages
        'org-babel-load-languages '((python . t) (R . t)))

Update: My next post shows how to call code in written in one language from code written in another language.

Related posts

Programming language popularity

Posted on by John

Here are two ways of measuring programming language popularity:

  1. Rank by number of questions tagged with that language on Stack Overflow
  2. Rank by number of project on GitHub using that language

According to this article, these two measures are well correlated.

I’d be skeptical of either metric by itself. A large number of questions on a language could indicate that it’s poorly documented, for example, rather than popular. And GitHub projects may not representative. But the two measures give similar pictures of the programming language landscape, so together they have more credibility. On the other hand, both measures are probably biased in favor newer languages.

The RedMonk Programming Language Rankings: February 2012

Example of not inverting a matrix: optimization

Posted on by John

People are invariably surprised when they hear it’s hardly ever necessary to invert a matrix. It’s very often necessary solve linear systems of the form Ax = b, but in practice you almost never do this by inverting A. This post will give an example of avoiding matrix inversion. I will explain how the Newton-Conjugate Gradient method works, implemented in SciPy by the function fmin_ncg.

If a matrix A is large and sparse, it may be possible to solve Ax = b but impossible to even store the matrix A-1 because there isn’t enough memory to hold it. Sometimes it’s sufficient to be able to form matrix-vector products Ax. Notice that this doesn’t mean you have to store the matrix A; you have to produce the product Ax as if you had stored the matrix A and multiplied it by x.

Very often there are physical reasons why the matrix A is sparse, i.e. most of its entries are zero and there is an exploitable pattern to the non-zero entries. There may be plenty of memory to store the non-zero elements of A, even though there would not be enough memory to store the entire matrix. Also, it may be possible to compute Ax much faster than it would be if you were to march along the full matrix, multiplying and adding a lot of zeros.

Iterative methods of solving Ax = b, such as the conjugate gradient method, create a sequence of approximations that converge (in theory) to the exact solution. These methods require forming products Ax and updating x as a result. These methods might be very useful for a couple reasons.

  1. You only have to form products of a sparse matrix and a vector.
  2. If don’t need a very accurate solution, you may be able to stop very early.

In Newton’s optimization method, you have to solve a linear system in order to find a search direction. In practice this system is often large and sparse. The ultimate goal of Newton’s method is to minimize a function, not to find perfect search directions. So you can save time by finding only approximately solutions to the problem of finding search directions. Maybe an exact solution would in theory take 100,000 iterations, but you can stop after only 10 iterations! This is the idea behind the Newton-Conjugate Gradient optimization method.

The function scipy.optimize.fmin_ncg can take as an argument a function fhess that computes the Hessian matrix H of the objective function. But more importantly, it lets you provide instead a function fhess_p that computes the product of the H with a vector. You don’t have to supply the actual Hessian matrix because the fmin_ncg method doesn’t need it. It only needs a way to compute matrix-vector products Hx to find approximate Newton search directions.

For more information, see the SciPy documentation for fmin_ncg.

More: Applied linear algebra