Thursday, March 31, 2016

What I was Reading Last Week

1. Evidence-based medicine has been hijacked: a report to David Sackett (John Ioannidis)

2. His interview at Retraction Watch

3. An a less pessimistic viewpoint (Steve Novella)

4. Another plea for Science-Based Medicine (David Gorski)

5. Massimi Pigliucci gives a Stoicism 101 talk (video)

6. What is the difference between Machine Learning and Data Science? (medium)

7. Family Tree versus phylogeny (nice cartoon inside)

8. Why Bill Nye was wrong when he dismissed Philosophy (patheos)

9. Andy Grove and the iPhone SE (stratechery)

10. A visual explanation of Markov Chains (Victor Powell)

Monday, March 28, 2016

Matrix Decompositions

I recently read G. W. Stewart's perspective "Decompositional Approach to Matrix Computations" published in Computing in Science in Engineering (2000).

He extols the virtues of this approach:
  • matrix decompositions solve many problems;
  • they are expensive to compute the first time, but can often be reused;
  • it facilitates rounding error analysis;
  • it fosters the development of general-purpose linear algebra software packages.
He goes on to describe six of the most famous matrix decompositions in his paper. Let me focus on the history of three of those six. It is interesting that much of the work on decomposition of matrices predates the relatively recent emergence of matrix algebra, which was fostered by digital computers.

1. Cholesky

For a positive definite matrix, \[A = LL^T.\] The legendary German mathematician Carl Friedrich Gauss (1777-1855) sketched out a method in 1809-1810.  It was more than  a century later that Andre-Louis Cholesky's variant was published posthumously by Benoit in 1924. Interestingly, Cholesky was a French military officer and mathematician, who used the decomposition in some of his surveying work, and Benoit was a fellow officer.

The closely related generalization, pivoted LU decomposition (\(PA = LU\)) has a more obscure history, although it does appear in Jacobi's work by 1857. Interestingly, pathological cases can be designed for which Gauss elimination with partial pivoting (the most popular strategy for LU decomposition) is unstable. However, the method works with remarkable versatility for most real-life problems, although why it is so successful is not completely understood.

2. QR decomposition

The QR decomposition, \[A = QR,\] first arose in early 1900 work on integral equations by Erhard Schmidt, a Baltic-German mathematician of the "Gram-Schmidt" orthogonalization fame. Alston Householder's method to triangularize a matrix using his namesake transformations was published in 1958, while Givens published  "Computation of plane unitary rotations transforming a general matrix to triangular form", also in 1958.

Both of them were at Oak Ridge National Laboratory around the time they made their seminal contributions. Incidentally, not only did their work on QR decomposition appear in same year (1958), they both also died in the same year (1993) 35 years later.

3. SVD

The singular value decomposition \[A = U \Sigma V^{T},\] can be thought of as a generalization of the spectral decomposition of symmetric matrices. Apparently, the SVD has a longer traceable history than the QR decomposition; it was introduced independently by Beltrami and Jordan in the mid 1870s. As Stewart says,
Together, Beltrami are Jordan are the progenitors of the SVD, Beltrami by virtue of first publication, and Jordan by the completeness and elegance of his treatment.
In 1965, Golub and Kahan published the reduction to the bidiagonal form. In a subsequent paper in 1970, Golub and Reinsch presented the (currently) most popular computational algorithm to perform SVD.

Wednesday, March 23, 2016

sed Converter for LaTeX Math Modes

In jupyter/iPython LaTeX markdown, I have many files that use the "$" and "$$" tags for math mode.

Sometimes, I liked to cut and paste stuff from there into Blogger. I have set my Blogger MathJax preferences to "slash-parentheses" and "slash-square-brackets" for "$" and "$$", respectively.

The following sed script lets me quickly do the conversion:

As as example consider the following:

Monday, March 21, 2016

QuickTip: Numpy logical AND

Consider the following logical operations in Matlab or Octave

> a=[1, 2, 3]
> b=[2, 2, 2]
> a >= 2 & b<= 2
 0   1   1

which might be what you expect. The conditions are applied element by element, and the output is a boolean array the same size as the input vectors.

In numpy:

>>> a=np.array([1,2,3])
>>> b=np.array([2,2,2])
>>> a >= 2 & b <= 2  
Traceback (most recent call last):
  File "", line 1, in
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

However, to get the expected response (array of bools) stick the two parts in parenthesis:

>>> (a >= 2) & (b <= 2)
array([False,  True,  True], dtype=bool)

Note the usual "and" operator doesn't work:

>>> (a >= 2) and (b <= 2)
Traceback (most recent call last):
  File "", line 1, in
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Saturday, March 19, 2016

What I Was Reading Last Week

1. Making inedible edible: GMO over the last 10,000 years (skepticalraptor)

2. Note to GMO skeptics: Get the freaking Science right! (Plato's footnote)

3. Maryam Namazie on Islamism (the Guardian)

4. On Sam Harris, Islamism, Free Speech and Hypocrisy (

5. Free speech is threatened on US college campuses (Intelligence Squared US)

6.  Another famed psychological truth falls? (Slate)

7. The menace of reproducibility (qz)

8. The irreproducibility of irreproducibility (Harvard Gazette)

9. Cancer: Fighting from the gut (Bloomberg)

10. Should academic literature be free? (NYT)

ReSpect: Program to extract Relaxation Time Spectra

During the last Society of Rheology meeting, I presented our work on extracting relaxation spectra from complex moduli. The goal is to extract the continuous relaxation spectrum \(h(\tau)\), which is given by :
\[\begin{align}G^{\prime}(\omega) & = \int\limits_{-\infty}^{\infty} \frac{\omega^2 \tau^2}{1 + \omega^2 \tau^2} \color{blue}{h(\tau)} d \ln \tau, \nonumber \\
G^{\prime\prime}(\omega) & = \int\limits_{-\infty}^{\infty} \frac{\omega \tau}{1 + \omega^2 \tau^2} \color{blue}{h(\tau)} d \ln \tau.
\end{align}\]Here are the slides from the presentation. The original version of this work appeared in the journal Applied Rheology, nearly two years ago.

Takeh A, Shanbhag S: A Computer Program to Extract the Continuous and Discrete Relaxation Spectra from Dynamic Viscoelastic Measurements, Appl. Rheol. 23 (2013) 24628.

The goal was to create a program to infer spectra with the design following attributes:
  • platform independence
  • ease-of-use "for an experimentalist"
  • transparent algorithm and implementation
  • free availability
  • continuous and discrete relaxation spectra
  • efficiency (~10 seconds to compute)
  • readability and extensibility of the code(*)
  • integrated graphics(*)
  • extension to modern multicore machines(*)
where * = optional.

Based on feedback I received in the months leading up to, and after the presentation, I put in a fair amount of work to improve the program, including bug-fixes, and better default parameter choices. I also added the capability of inferring the spectra from \(G(t)\) data as well. The program is hosted here.

My hope is to demonstrate a few cases here on the blog. 

Monday, March 14, 2016

How to Learn: GPS verus Map mentality.

A GPS is extremely convenient:
  • you can just get in the car, without any prior preparation
  • you get turn-by-turn instructions
  • if you lose your way, it recalculates an optimal path
  • there is no need to shuffle through print outs or maps, when you feel lost
One could approach a new class or subject with a "modular" GPS mentality. Simply focus on the task at hand. Don't think too much about what roads you took in the past, or need to take in the future. Let the instructor (or the GPS) take charge of all that.

Your job as the student (or driver) is to follow instructions as best as you can. If you follow them correctly, you will reach your destination, which in this context might involve earning a good grade. If you make a mistake, you can rely on the GPS or teacher to help you get back on.

While convenient and effective (if the only goal is a good grade), such an approach misses insights that might be available with a more struggle-prone Map mentality.

When you approach a new subject with a Map mentality, you expend effort to get a feel for the lay of the land; how are the important highways oriented? Are you headed north, south, east or west?

There is more contextual learning and reflection; an awareness of how newly learned material fits in with all that you have learned before. New knowledge is weaved into the fabric of prior concepts - leading to a richer latticework of connected ideas.

Of course, this dichotomy between GPS and Map mentalities is completely contrived -- there is no reason why you can't have them both. However, the convenience and apparent sufficiency of the GPS approach, often discourages a time-constrained student from looking at the Map or the big-picture.

In my opinion, such big-picture awareness is crucial for long-term education, as embodied in a quote I am fond of: "education is what remains, after all the details learned are forgotten!"

Thursday, March 10, 2016

Python > Matlab?

I find myself substituting Python (+Numpy+SciPy+SymPy) for applications, where I would have whipped out some quick Matlab/Octave code in the past. By no means do I consider myself a skilled python programmer, but almost everything I used to do in Octave/Matlab, I can do in python; and very often I can do it better.

As an analogy, Matlab is like an applicant who has a very high score on the quantitative section of the SAT; while Python is like an applicant who is nearly as good on quant, but also has excellent verbal, and writing skills.

Here are a set of articles that echo a similar sentiment:

1. Why use Python for Scientific Computing (Cyrille Rossant)
  • free and open-source
  • better and complete language
  • rich library set
  • plays nice with other languages and OS
2. Keep Calm and Code in Python (Lorena Barba)
  • great teaching language
  • very easy to get useful work done
  • good string manipulation
  • demand for python programmers
There are a lot of other links to python evangelists in the links above.

Personally, what I really like about python is:
  • the arbitrary constraint of requiring files and global functions to be "atomic"
  • I find visualization using matplotlib (with seaborn or ggplot) to be better and more configurable
  • same with 3D visualization using MayaVi2 for instance or arbitrary graphics
  • I just love jupyter notebooks. They help me keep code with documentation (which might include LaTeX markup). It has become my default mode of preparing new lectures.
  • python clearly has momentum

What are possible downsides of moving from Matlab to Python?
  • if you've used Matlab for a long time, you have to unlearn a few things
  • installation, while much simpler these days, can be intimidating for some
  • the python world currently is stuck between v2.x and v3.x
  • for a subset of tasks (especially linear algebra) I find Matlab more intuitive. For instance, 1d arrays are either column or row vectors; in numpy the distinction is fuzzy.
  • operator overloading: multiply matrix and vector in Matlab is just A*x; in python you have to say,x)

Friday, March 4, 2016


1. What would happen if GMO crops were banned (Steve Novella)
The authors of the current study demonstrate one key advantage of using GM traits – increased yield which results in lower land use for agriculture. They also point out that banning GMO will result in higher food prices, of $14-$24 billion per year in the US alone. Higher food prices disproportionately affect those with lower income.
2. A brutal takedown of "the Math Myth" (Keith Devlin)
The tragedy of The Math Myth is that Hacker is actually arguing for exactly the kind of life-relevant mathematics education that I and many of my colleagues have been arguing for all our careers. [...] Unfortunately, and I suspect because Hacker himself did not have the benefit of a good math education, his understanding of mathematics is so far off base, he does not recognize that the examples he holds up as illustrations of bad education only seem so to him, because he misunderstands them.
 3. John D. Cook excerpts an interesting take on rigor versus vigor in mathematics
I just started reading Frequency Analysis, Modulation and Noise by Stanford Goldman. The writing is strikingly elegant and clear. Here is a paragraph from the introduction.
Rigorous mathematics has a rightful place of honor in human thought. However, it has wisely been said that vigor is more important than rigor in the use of mathematics by the average man. In the particular case of this volume, the amount of rigor will be used that is necessary for a thorough understanding of the subject at hand by a radio engineer; but when it appears that rigor will confuse rather than clarify the subject for an engineer, we shall trust in the correctness of the results established by rigorous methods by the pure mathematicians and use them without the background of a rigorous proof.
In the class I currently teach, we try to adopt a similar posture. I also remembered a similar line is Sanjoy Mahajan's book " “Too much rigor teaches rigor mortis: the fear of making an unjustified leap even when it lands on the correct result.”

QuickTip: Save NumPy arrays to File as Column Vectors


import numpy as np

x=np.array([1, 2, 3])

print np.c_[x,y]
array([[1, 1],
       [2, 4],
       [3, 9]])

So if you want to save two arrays (of the same size) to file then,

np.savetxt('xy.dat', np.c_[x,y], fmt='%3.4f')