Saturday, March 28, 2015

For the love of solitude

Last week I spent some time in a small cabin ("chalet") in the woods. At the time it was completely deserted: the silence was delicious. Finally I could breath. Finally my thoughts were my own. I am always struck in these moments, first, how essential they are to someone of my temperament, and second, how it clears the mind so that the real thinking may finally begin.

Such solitude is increasingly hard to find. I point out in another article how the world is becoming a panopticon. Soon there will be satellites with sufficient resolution and coverage that they can observe us in real time. Not even walls will be enough to conceal us: the combination of penetrating microwaves, tomography and synthetic apertures will soon allow us (or rather our nosy leaders) to map the insides of buildings in 3-D from space.

I think I finally understand what's going on here; why I am driven to find more and more extreme isolation, so much so that loneliness and anxiety frequently overcomes any peace-of-mind that might be gained from the endeavour.

When you observe something, it changes.  Many interpretations of quantum mechanics have been attempted, many of them rather flaky: something is not there when you are not looking at it; there is a mystic union between the observed and the observer such that they cannot be distinguished.  Lets not even get into all the different "many-worlds" hypotheses.

No, it is much simpler than that, and from the mathematics, undeniable.  Chances are you can walk into any toy store today and buy a simple device consisting of an array of needles free to slide within a series of holes.  Using this device, you can take a temporary cast of your face (or any other object for that matter) by simply pressing into it.  Now at the same time your face is pushing these needles outwards, the needles are creating pock-mark depressions on your skin.  Granted, because skin is elastic, it will almost certainly spring back to it previous form, although you might feel it for some time afterwards.

Rays of light are just like those needles.  When you take leave to walk in the woods or cycle in the mountains, you are reclaiming your thoughts.  You are becoming fully yourself again, because as long as others watch you, your thoughts are not your own.

It might well be my epitaph.  Such solitude is no longer supportable.  I guess I consider myself a dying breed.

How to deal with objects

In my coding standards, I state that I strive for a predominantly functional programming model, even when using an imperative, procedural language such as C++. All functions as well as main routines should, in the ideal case, take a set of inputs and return a set of outputs while avoiding any kind of side effect. This leaves open the question of how to deal with objects which I use quite often. Object-oriented programming doesn't quite fit the functional paradigm, so here is my attempt to lay out an approach consistent with my earlier standards.

Class data ("fields") can be divided into two, broad categories: those that define the object and workspace variables. As an example of the first type, consider a matrix class. A matrix object would be defined by its dimensions and the values of each individual element. These can be modified by the class methods such that the method can be considered an operator on the object. So for instance the matrix class might have a "transpose" or "invert" method. More basic operators might include subscript assignment, permuting rows or columns, and in particular, reading from and writing to a file.

 The other type of class data are workspace variables. Using the matrix example: internally the inverse method might comprise two steps--the first to decompose the matrix and the second to perform some other operation on the decomposition. The decomposition is stored inside of the object. This approach gives a minor performance advantage: if we know in advance the size of the workspace variable, it can be pre-allocated as soon as the the object has been defined. It also simplifies memory management.

 My current view on object use and modification is rather extreme: data fields should be filled completely at or near the moment of creation and from that moment on the object should be considered more-or-less immutable. Methods may take other objects as arguments and return newly-created objects, but modifications to an existing object should be avoided in favour of deleting it and creating a new one in its place.

Workspace variables in particular should be avoided--they are, in fact, global variables in disguise. This is especially so if they are being used in a context that doesn't modify the more essential object data and is meant to keep the re-entrant property and make the code thread-safe. Operators that modify the object, especially if they are simple such as subscript assignment, can always be made atomic. By contrast, workspace variables by their very nature frequently span multiple methods. Computers are now so fast and memory management so efficient that the prior justification of improving efficiency is nullified. Enabling the program to run on multiple processors will dwarf any marginal gains that could be had with the use of workspace variables.

I admit that most of my current code isn't anywhere close to conforming to these standards although I rarely use workspace variables and have started to actively remove them from my O-O code.

Tuesday, February 17, 2015


A quine is a program that prints out it's own listing. In other words, if a it's a C program, you can do something like this:

$ gcc quine.c
$ ./a.out > quine2.c
$ gcc
$ ./a.out > quine3.c
$ diff quine3.c quine.c

Writing a quine is not as trivial as it sounds.  Here is the quine I wrote in C which takes advantage of the format codes in the printf statement to insert a partial listing of the program into itself:

#include <stdio.h> int main() {char listing[200]="#include <stdio.h> %cint main() {char listing[200]=%c%s%c; char delim=10; char delim2=34; printf(listing, delim, delim2, listing, delim2);}"; char delim=10; char delim2=34; printf(listing, delim, delim2, listing, delim2);}

Sunday, January 18, 2015

The Overseer

One of the more interesting things I did this summer was visit an old-growth forest.  While the scale of these trees was truly astounding, what struck me most was just how little of the tree was canopy and how much was trunk.  Most of the mass of these trees is there to hold up other parts of the tree!  Only a tiny fraction is actually devoted to gathering sunlight.

Moreover, there will be an equilibrium density for the trees.  If they are clustered together too tightly, then one or more of their number will fall.  If they are too sparse, new trees will grow up in the gaps.  In other words, in their natural state, they live a marginal existence: clinging to the edge of life.

That is, unless there is someone tending the forest who culls the trees as necessary.  Then those remaining can grow up healthier and stronger, with thicker trunks relative to their height and a larger, broader canopy.

Monday, January 12, 2015

What is spirit?

In a two or more posts, I've talked a little bit about "spirit".  If you're a skeptic and you've been reading this, you probably think, that's nonsense, there's no such thing as spirit, it doesn't exist.  Probably not, but then again, there isn't any such thing as energy either, at least not in the sense that you can measure it directly.  It is a derived quantity: it always has to be calculated from two or more other quantities.  Nonetheless it's still an extremely useful quantity in physics.

In particular, I wanted to clarify some statements I made in the post called, My Dream. I thought about simply writing an addendum, but why not create a brand new post?  In that post I talk about negative spirit flowing out of Washington D.C. through a network of trails leading, ultimately, to long-distance scenic trails like the Appalachian Trail.  How does that work?  What is the mechanism behind this?  Simple.  The people walking (or cycling) towards the city are smiling.  Those moving away are not.

I experienced this directly.  I picked up two hikers from the A.T. and let them sleep over in my apartment.  The next day, we made arrangements to meet after I had finished work and they had finished sight-seeing.  Later, I began to worry that I would having trouble finding them since we had agreed to meet at one of the busiest subway stations during rush hour.  It turns out they were easy to spot: they were the only ones smiling.

Sunday, January 11, 2015

Are miracles possible?

What about the double-slit experiment?  That's kind of a miracle isn't it?  Point an electron gun at a pair of slits and then look at where they land on a screen behind.  It doesn't matter how fast you fire the electrons, they will, over time, form an interference pattern.  In other words, they are behaving like waves.  OK, so now you're wondering, if they are behaving like waves, yet they're being fired one at a time, what is interfering with what?  Does each electron interfere with itself, or do they interfere with each other, somehow "seeing" the other electrons both forward and backward in time?

Lets design an experiment to test these ideas.  Now you flash a strobe light every time you fire an electron so you can figure out through which slit it passes.  The interference pattern disappears.

I've heard any number of interpretations of quantum mechanics: An object isn't really there when you aren't observing it.  There is a mystical connection between observer and observed.  The distinction between subject and object is an illusion. And so on.

But at least in this case, I think the interpretation is much, much simpler.  When you observe something, it changes.  Suppose instead of using the term observe, we substitute the word probe.  To probe is almost by definition to interfere with since a probe (as a noun) is almost always something material and solid.  Granted, if you drop a sounder into the ocean to measure depth, for instance, the ocean floor won't change much in response, certainly less than measurement error or the local variance.  But if you fire a photon at a sub-atomic particle, it will change the velocity of that particle, sometimes by a large amount.

The act of observation is the act of touching.  When something is touched, it will always move or change in response.

The analogy to investigations of so-called "psi" phenomenon is (or should be) obvious.  Lets use prayer as an example.  I have at various times in my life, prayed with considerable regularity.  While I no longer identify as Christian, I still pray, though less often than I used to.  It was a great sense of relief when I stopped attending church because I could finally admit to myself that I am completely ignorant about spiritual matters.  But then again, the whole point of a Supreme Being is that He or She or It is completely beyond us.

I don't know how many times I've heard of studies purporting to demonstrate the effectiveness of prayer.  I don't know how many times I've heard skeptics try to debunk such claims and reference their own studies.  But after all, if spiritual configurations (for lack of a better word) are as delicate as sub-atomic particles, should we be surprised if close observation destroys them?  Should we be surprised if they behave in the same way as the double slit experiment and lose their magic as soon as they are seriously scrutinized?

I don't know the answer to these questions, because there is no answer.  I guess it just comes down to faith.

Saturday, January 10, 2015

Fun with matrix operators

Now that I have more time, I can turn over ideas more carefully.  I'm not so focused on results. Take matrix multiplication.  It is an associative operator:


It's a property I use all the time, but typically never consciously think about. For instance, in my sparse matrix array classes, the sparse matrices are stored in reverse order because they are multiplied in reverse order.  This is for efficiency reasons.  If you have a set of matrices, {Ai}, and you want to multiply them through and then multiply the result with a vector, as follows:

A1A2A3 ... An-1An v                                (1)

then this:

A1(A2(A3( ... (An-1(An v)) ... )))             (2)

is more efficient than multiplying them through in the order in which they are written:

(( ... ((A1A2)A3) ... An-1)An v)               (3)

because each intermediate result is a vector, rather than a matrix.  Similarly, if v is a full matrix instead of a vector, (2) is still more efficient than (3) because multiplying a sparse matrix with a full one is more efficient than multiplying two sparse matrices.

Is associativity related to linearity?  No, because we can define a new, linear operator, lets call it, #:

A#B = ABT = Σk aikbik

that loses the associative property.  Now lets rewrite (1) in terms of #:

A1#(A2#(A3( ... (An-1#(An vT)T)T ... )T)T)T

                     A1#vT#An#An-1# ... #A3#A2

                    = (vT#An#An-1# ... #A1#A2#A1)T                  (4)


(A#B)T = B#A

So why is this useful?  From a strictly mathematical standpoint, the new operator is quite useless: it loses the useful and important property of associativity, without producing any new advantages.  Numerically, however, it is extremely useful: multiplying a sparse matrix with the transpose of a sparse matrix is much more straightforward than straight multiplication since you move through the non-zero elements of both matrices in order.  This is assuming that both are sorted by subscript in row-major order.

Suppose, once again, that v is a matrix.  Even if v is sparse, it's usually more efficient to convert it to a full matrix and multiply through using (2) since the result will become full after only a few iterations.  Not always: if all the matrices are tri-diagonal, for instance, then the result will also be tri-diagonal.

So lets assume that not only is v a sparse matrix, but that the result will also be sparse.  Clearly, (4) is now more efficient than (2) since we now need only 2 matrix transpositions instead of n, where n is the number of A's. Since matrix multiplication is O(n2) (where n is the number of non-zero elements in each matrix, assuming it is about the same for each) while matrix transposition is O(n log n) (for a generalized sparse matrix: swap the subscripts, then re-sort them in row-major order) the gains will be quite modest.

When constructing parallel algorithms, it might be better to stick with (2) since we can multiply through in just about any order we want: e.g., multiply, in order, every second pair of consecutive matrices, then take the results and again multiply by consecutive pairs and so on.

I mention multiplying an array of sparse matrices with a full matrix.  This is also an interesting case, although mostly from a numerical perspective.  The intuitive approach would be to have two full matrices: one for the result, and one for the multiplicand. We multiply each of the A's in turn, swapping the result with the multiplicand in between.  A second approach would be to decompose the multiplicand into columns and multiply each of the columns in turn, all the way through.

It turns out that for large problems, the latter approach is much, much faster, even though it violates at least one optimization principle. If we have two loops, both of which have the same limits, convert them to one loop.  For example:

for (int i=0; i < n; i++) a[i]=b[i]*b[i];
for (int i=0; i < n; i++) r+=a[i];

should be converted to:

for (int i=0; i < n; i++) {

The column-by-column method is better because it uses less memory and because it operates on a smaller chunk of memory for longer.  Whereas the first method accesses a whole full matrix at each multiplication step, the second method accesses only one column of the matrix. Therefore there is less paging, either between virtual memory and RAM or between core memory and cache.  Like multiplication between two sparse matrices, this method is also more efficient when multiplying with the transpose, assuming row-major storage order once again for both the sparse matrices and the full matrices.

It does well to remind us that while high-level languages like C and C++ do a good job isolating the programmer from the hardware, you don't always end up with the best code by considering the computer as this perfect, abstract machine, instead of a real one with real limitations and real, hacked-together solutions and work-arounds.  

It is worth comparing the code for both methods, so I will post it in full.  Here is method #1:

void sparse_array::mat_mult(real **cand,

                real **result) {
  real **med;
  real **swap;

  med=allocate_matrix(m, m);
  sparse_a[0]->mat_mult(cand, result, m);
  for (long i=1; i
    sparse_a[i]->mat_mult(result, med, m);

  if (nsparse % 2 == 0) {
    copy_matrix(result, med, m, m);

Here is method #2, applied, however, to the transpose of the multiplicand:

void sparse_array::mat_mult_tt(real **cand,
                real **result) {
  printf("%6.1f%%", 0.);
  for (index_t i=0; i
    printf("\b\b\b\b\b\b\b%6.1f%%", 100.*i/m);
    vect_mult(cand[i], result[i]);
  printf("\b\b\b\b\b\b\b%6.1f%%\n", 100.);

There are a few of things worth noting.  In the first method, if the number of sparse matrices is odd, the result and the multiplicand are swapped an odd number of times.  For efficiency, only the memory locations are swapped, not the actual values so the memory locations for the multiplicand and the result end up swapped.  Therefore we need to copy the values back to the correct location in memory before the function can return.  In the second method, we print out a crude progress meter in the form of the percentage of rows that have been multiplied through.

The second version is not only much simpler, it is also easier to parallelize since each row multiplication is completely independent.