Walking Randomly

Last week I gave a live demo of the IPython notebook to a group of numerical analysts and one of the computations we attempted to do was to solve the following linear system using Numpy’s solve command.

Now, the matrix shown above is singular and so we expect that we might have problems. Before looking at how Numpy deals with this computation, lets take a look at what happens if you ask MATLAB to do it

>> A=[1 2 3;4 5 6;7 8 9];
>> b=[15;15;15];
>> x=A\b
Warning: Matrix is close to singular
or badly scaled. Results may be
inaccurate. RCOND =  1.541976e-18. 
x =
  -39.0000
   63.0000
  -24.0000

MATLAB gives us a warning that the input matrix is close to being singular (note that it didn’t actually recognize that it is singular) along with an estimate of the reciprocal of the condition number. It tells us that the results may be inaccurate and we’d do well to check. So, lets check:

>> A*x

ans =
   15.0000
   15.0000
   15.0000

>> norm(A*x-b)

ans =
2.8422e-14

We seem to have dodged the bullet since, despite the singular nature of our matrix, MATLAB has able to find a valid solution. MATLAB was right to have warned us though…in other cases we might not have been so lucky.

Let’s see how Numpy deals with this using the IPython notebook:

In [1]:
import numpy
from numpy import array
from numpy.linalg import solve
A=array([[1,2,3],[4,5,6],[7,8,9]])
b=array([15,15,15])
solve(A,b)

Out[1]:

array([-39.,  63., -24.])

It gave the same result as MATLAB [See note 1], presumably because it’s using the exact same LAPACK routine, but there was no warning of the singular nature of the matrix.  During my demo, it was generally felt by everyone in the room that a warning should have been given, particularly when working in an interactive setting.

If you look at the documentation for Numpy’s solve command you’ll see that it is supposed to throw an exception when the matrix is singular but it clearly didn’t do so here. The exception is sometimes thrown though:

In [4]:

C=array([[1,1,1],[1,1,1],[1,1,1]])
x=solve(C,b)

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
 in ()
      1 C=array([[1,1,1],[1,1,1],[1,1,1]])
----> 2 x=solve(C,b)

C:\Python32\lib\site-packages\numpy\linalg\linalg.py in solve(a, b)
    326     results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
    327     if results['info'] > 0:
--> 328         raise LinAlgError('Singular matrix')
    329     if one_eq:
    330         return wrap(b.ravel().astype(result_t))

LinAlgError: Singular matrix

It seems that Numpy is somehow checking for exact singularity but this will rarely be detected due to rounding errors. Those I’ve spoken to consider that MATLAB’s approach of estimating the condition number and warning when that is high would be better behavior since it alerts the user to the fact that the matrix is badly conditioned.

Thanks to Nick Higham and David Silvester for useful discussions regarding this post.

Notes
[1] – The results really are identical which you can see by rerunning the calculation after evaluating format long in MATLAB and numpy.set_printoptions(precision=15) in Python

Along with a colleague, I’ve been playing around with Anaconda Python recently and am very impressed with it. At the time of writing, it is at version 1.7 and comes with Python 2.7.5 by default but you can install Python 3.3 using their conda package manager.  After you’ve installed Anaconda, just start up a Windows command prompt (cmd.exe) and do

conda update conda
conda create -n py33 python=3.3 anaconda

It will chug along for a while, downloading and installing packages before leaving you with a Python 3.3 environment that is completely separated from the default 2.7.5 environment. All you have to do to activate Python 3.3 is issue the following command at the Windows command prompt

activate py33

To demonstrate that the standard anaconda build remains untouched, launch cmd.exe, type ipython and note that you are still using Python 2.7.5

Microsoft Windows [Version 6.1.7601]
Copyright (c) 2009 Microsoft Corporation.  All rights reserved.

C:\Users\testuser>ipython
Python 2.7.5 |Anaconda 1.7.0 (64-bit)| (default, Jul  1 2013, 12:37:52) [MSC v.1500 64 bit (AMD64)]
Type "copyright", "credits" or "license" for more information.

IPython 1.0.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]:

Exit out of ipython and activate the py33 environment before launching ipython again. This time, note that you are using Python 3.3.2

In [1]: exit()

C:\Users\testuser>activate py33
Activating environment "py33"...

[py33] C:\Users\testuser>ipython
Python 3.3.2 |Anaconda 1.7.0 (64-bit)| (default, May 17 2013, 11:32:27) [MSC v.1500 64 bit (AMD64)]
Type "copyright", "credits" or "license" for more information.

IPython 1.0.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]:

In a recent blog-post, John Cook, considered when series such as the following converged for a given complex number z

z₁ = sin(z)
z₂ = sin(sin(z))
z₃ = sin(sin(sin(z)))

John’s article discussed a theorem that answered the question for a few special cases and this got me thinking: What would the complete set of solutions look like? Since I was halfway through my commute to work and had nothing better to do, I thought I’d find out.

The following Mathematica code considers points in the square portion of the complex plane where both real and imaginary parts range from -8 to 8. If the sequence converges for a particular point, I colour it black.

LaunchKernels[4]; (*Set up for 4 core parallel compute*)
ParallelEvaluate[SetSystemOptions["CatchMachineUnderflow" -> False]];
convTest[z_, tol_, max_] := Module[{list},
  list = Quiet[
    NestWhileList[Sin[#] &, z, (Abs[#1 - #2] > tol &), 2, max]];
  If[
   Length[list] < max && NumericQ[list[[-1]]]
   , 1, 0]
  ]
step = 0.005;
extent = 8;
AbsoluteTiming[
 data = ParallelMap[convTest[#, 10*10^-4, 1000] &,
    Table[x + I y, {y, -extent, extent, step}, {x, -extent, extent,
      step}]
    , {2}];]
ArrayPlot[data]

I quickly emailed John to tell him of my discovery but on actually getting to work I discovered that the above fractal is actually very well known. There’s even a colour version on Wolfram’s MathWorld site.  Still, it was a fun discovery while it lasted

Other WalkingRandomly posts like this one:

• Complex PowerTowers
• Quadraflakes, Pentaflakes, Hexaflakes and more

Right now it’s packaging season (not the official term!) at my university–a time of year when IT staff have to battle with silent installers, SCCM, MSI creation and Adminstudio in order to create the student desktop image for the next academic year.  Packaging season makes me grouchy, it makes me work late and it makes me massively over-react to every minor installation issue caused by software vendors. Right now, however, I am not grouchy because of packaging season..I am grouchy because of concurrent network licensing (or floating licensing if you are Wikipedia).

I like network licenses…they make many aspects of my job easier but they the way they are implemented by some software vendors causes them to be a pain.  Over the years, I’ve bothered many a support-desk with my network license tails of woe and thought that I would collate them all together in one blog post.

The more of these things your software does, the more pain you cause me and my colleagues.

1. You don’t use standard FlexLM/FlexNet. 

Like it or loathe it, FlexLM is used by the vast majority of software vendors out there. We run license servers that host dozens of FlexLM based applications and we have got the administration of these down to a fine art.  In fact, we’ve replaced the vast majority of the process with a script. If an application uses FlexLM, system administration and license monitoring is bordering on the trivial for us now.  The further you stray away from FlexLM, the more difficult our job becomes.

One thing guaranteed to ruin my day is a call from a vendor I’ve worked with for years who says ‘Great news Mike, we’re ditching FlexLM for our own, in house license system.’  Fabulous! Rather than use our lovely, generic, one-size fits all scripts, we are going to have to do a load of extra work and testing just for you.  I look forward to all the new and interesting bugs your system will generate.

2. You don’t support redundant license servers.

The idea behind redundant license servers is this:  Instead of your application relying on just one machine, it relies on three; only two of which have to be operational at any one time.  This gives us resiliency and resiliency is a good thing when you are teaching a lab with 120 students in it.

I’ll keep this simple.  If you don’t support redundant license servers, it means that you don’t believe that your software is important. It tells me that you are just playing at being a big, grown up piece of software but you don’t really think anyone will take you seriously.

3. You support redundant license servers but can’t select them via the installer.

At install time, there is no option to give three severs. The user can only give one. You then expect the user to copy a pre-prepared license file that has details of all three servers as a post-installation step.

What usually happens here is that users give the primary license server, find that the application will launch and stop reading the installation instructions.  At some point in the future, we take down the primary license server for maintenance and the vast majority of self-serve installations break!

4. You use the LM_LICENSE_FILE environment variable

The problem is, so does everyone else. We end up with a situation where the LM_LICENSE_FILE variable is pointing at several license servers and some client applications really don’t like that. Be a good FlexLM citizen and use a vendor specific environment variable.  For example, MATLAB uses MLM_LICENSE_FILE and I could give them a big hug just for that!

5. You ‘helpfully’ tell the user when the license is about to expire.

I’ve moaned about this before. 1000 users panicking and emailing the helpdesk…lovely!  Bonus points are awarded if you don’t allow any supported way of switching these warnings off.

6. Your new license doesn’t support old clients.

This should speak for itself and happens more than I’d like.  We can’t upgrade the entire estate instantaneously and even if we could, we probably wouldn’t want to.  Some users, for one reason or another, cling to old versions of your software like grim death. They don’t care that there is a new shiny version available, all they know is that I broke their application and they hate me for it.

When we discover that old versions of your application will stop working, it also delays roll out of the new version since we have to do a lot of user-communications and ensure that nothing mission-critical will stop working.  This makes power-users of your application hate me because they want the new shiny version.

7. You don’t have a silent installer.

Strictly speaking not a network license moan but closely related.  We use network licensing because we deploy your software to lots of machines.  When I say ‘lots’ I mean thousands.  It turns out, however, that you don’t support scripted installation (sometimes called ‘silent installation’ or ‘unattended installation’).  This means that your software is a lot more difficult to deploy than your competitor!  I’m now a big fan of your competitor!

8. You have a silent installer but it’s a bad one.

If I install manually, via point and click, I can configure every aspect of your application.  Your silent installer, on the other hand, is just a /S switch that does a default install…no configuration possible.  Bonus points for ‘silent’ installers that include pop-up dialogue boxes that can’t be switched off.

While on the topic of silent installation, can I just ask that you directly support deployment by SCCM on Windows please?  It will help with next year’s packaging season big time!

Cheers,

Mike

The first stable version of KryPy was released in late July.  KryPy is “a Python  module for Krylov subspace methods for the solution of linear algebraic systems. This includes enhanced versions of CG, MINRES and GMRES as well as methods for the efficient solution of sequences of linear systems.”

Here’s a toy example taken from KryPy’s website that shows how easy it is to use.

from numpy import ones
from scipy.sparse import spdiags
from krypy.linsys import gmres

N = 10
A = spdiags(range(1,N+1), [0], N, N)
b = ones((N,1))

sol = gmres(A, b)
print (sol['relresvec'])

Thanks to KryPy’s author, André Gaul, for the news.

While working on someone’s MATLAB code today there came a point when it was necessary to generate a vector of powers.  For example, [a a^2 a^3….a^10000] where a=0.999

a=0.9999;
y=a.^(1:10000);

This isn’t the only way one could form such a vector and I was curious whether or not an alternative method might be faster. On my current machine we have:

>> tic;y=a.^(1:10000);toc
Elapsed time is 0.001526 seconds.
>> tic;y=a.^(1:10000);toc
Elapsed time is 0.001529 seconds.
>> tic;y=a.^(1:10000);toc
Elapsed time is 0.001716 seconds.

Let’s look at the last result in the vector y

>> y(end)
ans =
   0.367861046432970

So, 0.0015-ish seconds to beat.

>> tic;y1=cumprod(ones(1,10000)*a);toc
Elapsed time is 0.000075 seconds.
>> tic;y1=cumprod(ones(1,10000)*a);toc
Elapsed time is 0.000075 seconds.
>> tic;y1=cumprod(ones(1,10000)*a);toc
Elapsed time is 0.000075 seconds.

soundly beaten! More than a factor of 20 in fact. Let’s check out that last result

>> y1(end)
ans =
   0.367861046432969

Only a difference in the 15th decimal place–I’m happy with that. What I’m wondering now, however, is will my faster method ever cause me grief?

This is only an academic exercise since this is not exactly a hot spot in the code!

In common with many higher educational establishments, the University I work for has site licenses for a wide variety of scientific software applications such as Labview, MATLAB, Mathematica, NAG, Maple and Mathcad— a great boon to students and researcher who study and work there.   The computational education of a typical STEM undergraduate will include exposure to at least some of these systems along with traditional programming languages such as Java, C or Fortran and maybe even a little Excel when times are particularly bad!

Some argue that such exposure to multiple computational systems is a good thing while others argue that it leads to confusion and a ‘jack of all trades and master of none situation.’  Those who take the latter viewpoint tend to want to standardize on one system or other depending on personal preferences and expertise.

MATLAB, Python, Fortran and Mathematica are a few systems I’ve seen put forward over the years with the idea being that students will learn the basics of one system in their first few weeks and then the entire subject curriculum will be interwoven with these computational skills.  In this way, students can use their computational skills as an aid to deeper subject understanding without getting bogged down with the technical details of several different computational systems.  As you might expect, software vendors are extremely keen on this idea and will happily parachute in a few experts to help universities with curriculum development for free since this will possibly lead to their system being adopted.

Maybe we’ll end up with electrical engineers who’ve only ever seen Labview, mathematicians who’ve only ever used Maple, mechanical engineers who’ve only ever used MATLAB and economists who can only use Excel.  While the computational framework(s) used to teach these subjects are less important than the teaching of the subject itself, I firmly believe that part of a well-rounded, numerate education should include exposure to several computational systems and such software mono-cultures should be avoided at all costs.

Part of the reason for writing this post is to ask what you think so comments are very welcome.

I am currently in the market for a new thin and light laptop; the things that everyone seems to be calling ultrabooks these days.  Here’s what I’ve considered so far and the major good/bad points from my point of view.

• 2013 Macbook Air: Good: Fantastic hardware specs including Haswell and Intel HD 5000. I’ll get an educational discount. Bad: No touchscreen, I’m not interested in OS X. No tablet mode.
• Lenovo Yoga 11s.  Good: Both a tablet and a laptop but it’s a laptop first and a tablet second.  Bad: Only 3rd generation Intel.  Weird to hold in laptop mode because your fingers are mashing the keyboard.
• Haswell based Dell XPS 12.  Good:Haswell. Beautiful mechanism for switching between tablet/laptop. Bad: Only Intel HD 4400 graphics compared to the Air’s Intel HD 5000.  I’ve heard that the connection between computer and screen is unencrypted WiDi. Possibly just FUD but its putting me off right now.
• Samsung Ativ Book 9 Plus Good: Haswell. Amazing display resolution. Bad: Much higher resolution than the Macbook Air but with weaker GPU (Intel HD 4400)–that surely can’t be a good combination? No tablet mode. Win 8 doesn’t handle such high  resolution screens well but should be fixed in 8.1.
• Samsung Ativ Q: Good: Haswell. Amazing display resolution. Tablet mode.Win 8+Android.  Bad: No trackpad. Looks horrible in laptop mode–seems to be a tablet first and laptop second.  Like the Book 9 plus, it has a relatively weak GPU (Intel HD 4400) to drive all those pixels.

I think that at the moment, what I want is the form factor and weight of the Dell XPS12 with the resolution of the Samsung and the hardware specs of the Air.  I’m open to suggestions though.

A colleague recently sent me this issue. Consider the following integral

Attempting to evaluate this with Mathematica 9 gives 0:

 f[a_, b_] := Exp[I*(a*x^3 + b*x^2)];
Integrate[f[a, b], {x, -Infinity, Infinity}, 
 Assumptions -> {a > 0, b \[Element] Reals}]

Out[2]= 0

So, for all valid values of a and b, Mathematica tells us that the resulting integral will be 0.

However, Let’s set a=1/3 and b=0 in the function f and ask Mathematica to evaluate that integral

In[3]:= Integrate[f[1/3, 0], {x, -Infinity, Infinity}]

Out[3]= (2*Pi)/(3^(2/3)*Gamma[2/3])

This is definitely not zero as we can see by numerically evaluating the result:

In[5]:=
(2*Pi)/(3^(2/3)*Gamma[2/3])//N

Out[5]= 2.23071

Similarly if we put a=1/3 and b=1 we get another non-zero result

In[7]:= Integrate[f[1/3, 1], {x, -Infinity, Infinity}] // InputForm

Out[7]=2*E^((2*I)/3)*Pi*AiryAi[-1]

In[8]:=
2*E^((2*I)/3)*Pi*AiryAi[-1]//N

Out[8]= 2.64453 + 2.08083 I

We are faced with a situation where Mathematica is disagreeing with itself.  On one hand, it says that the integral is 0 for all a and b but on the other it gives non-zero results for certain combinations of a and b.  Which result do we trust?

One way to proceed would be to use the NIntegrate[] function for the two cases where a and b are explicitly defined.  NIntegrate[] uses completely different algorithms from Integrate.  In particular it uses numerical rather than symbolic methods (apart from some symbolic pre-processing).

NIntegrate[f[1/3, 0], {x, -Infinity, Infinity}]

gives 2.23071 + 0. I and

NIntegrate[f[1/3, 1], {x, -Infinity, Infinity}]

gives 2.64453 + 2.08083 I

What we’ve shown here is that evaluating these integrals using numerical methods gives the same result as evaluating using symbolic methods and numericalizing the result.  This gives me some confidence that its the general, symbolic evaluation that’s incorrect and hence I can file a bug report with Wolfram Research.

Maple 17.01 on the general problem

Since my University has just got a site license for Maple, I wondered what Maple 17.01 would make of the general integral. Using Maple’s 1D input/output we get:

> myint := int(exp(I*(a*x^3+b*x^2)), x = -infinity .. infinity);

myint := (1/12)*3^(1/2)*2^(1/3)*((4/3)*Pi^(5/2)*(((4/27)*I)*b^3/a^2)^(1/3)*exp(((2/27)*I)*b^3/a^2)
*BesselI(-1/3, ((2/27)*I)*b^3/a^2)+((2/3)*I)*Pi^(3/2)*3^(1/2)*b*2^(2/3)
*hypergeom([1/2, 1], [2/3, 4/3],((4/27)*I)*b^3/a^2)/(-I*a)^(2/3)-(8/27)*2^(1/3)*Pi^(5/2)*b^2*
exp(((2/27)*I)*b^3/a^2)*BesselI(1/3, ((2/27)*I)*b^3/a^2)/((-I*a)^(4/3)*(((4/27)*I)*b^3/a^2)^(1/3)))
/(Pi^(3/2)*(-I*a)^(1/3))+(1/12)*3^(1/2)*2^(1/3)*((4/3)*Pi^(5/2)*(((4/27)*I)*b^3/a^2)
^(1/3)*exp(((2/27)*I)*b^3/a^2)*BesselI(-1/3, ((2/27)*I)*b^3/a^2)+((2/3)*I)*Pi^(3/2)*3^(1/2)*b*2^(2/3)
*hypergeom([1/2, 1], [2/3, 4/3], ((4/27)*I)*b^3/a^2)/(I*a)^(2/3)-(8/27)*2^(1/3)*Pi^(5/2)*b^2*
exp(((2/27)*I)*b^3/a^2)*BesselI(1/3, ((2/27)*I)*b^3/a^2)/((I*a)^(4/3)*(((4/27)*I)*b^3/a^2)^(1/3)))
/(Pi^(3/2)*(I*a)^(1/3))

That looks scary! To try to determine if it’s possibly a correct general result, let’s turn this expression into a function and evaluate it for values of a and b we already know the answer to.

>f := unapply(%, a, b):
>result1:=simplify(f(1/3,1));

result1 := (2/27)*3^(1/2)*Pi*exp((2/3)*I)*((-(1/3)*I)^(1/3)*3^(2/3)*(-1)^(1/6)*BesselI(-1/3, (2/3)*I)
+3*(-(1/3)*I)^(2/3)*BesselI(-1/3, (2/3)*I)+3*(-(1/3)*I)^(2/3)*BesselI(1/3, (2/3)*I)-BesselI(1/3, (2/3)*I)
*3^(1/3)*(-1)^(1/3))/(-(1/3)*I)^(2/3)

evalf(result1);
2.644532853+2.080831872*I

Recall that Mathematica returned 2*E^((2*I)/3)*Pi*AiryAi[-1]=2.64453 + 2.08083 I for this case. The numerical results agree to the default precision reported by both packages so I am reasonably confident that Maple’s general solution is correct.

Not so simple simplification?

I am also confident that Maple’s expression involving Bessel functions is equivalent to Mathematica’s expression involving the AiryAi function. I haven’t figured out, however, how to get Maple to automatically reduce the Bessel functions down to AiryAi. I can attempt to go the other way though. In Maple:

>convert(2*exp((2*I)/3)*Pi*AiryAi(-1),Bessel);

(2/3)*exp((2/3)*I)*Pi*(-1)^(1/6)*(-BesselI(1/3, (2/3)*I)*(-1)^(2/3)+BesselI(-1/3, (2/3)*I))

This is more compact than the Bessel function result I get from Maple’s simplify so I guess that Maple’s simplify function could have done a little better here.

Not so general general solution?

Maple’s solution of the general problem should be good for all a and b right?  Let’s try it with a=1/3, b=0

f(1/3,0);
Error, (in BesselI) numeric exception: division by zero

Uh-Oh! So it’s not valid for b=0 then! We know from Mathematica, however, that the solution for this particular case is (2*Pi)/(3^(2/3)*Gamma[2/3])=2.23071. Indeed, if we solve this integral directly in Maple, it agrees with Mathematica

>myint := int(exp(I*(1/3*x^3+0*x^2)), x = -infinity .. infinity);

myint := (2/9)*3^(5/6)*(-1)^(1/6)*Pi/GAMMA(2/3)-(2/9)*3^(5/6)*(-1)^(5/6)*Pi/GAMMA(2/3)

>simplify(myint);
(2/3)*3^(1/3)*Pi/GAMMA(2/3)

>evalf(myint);
2.230707052+0.*I

Going to back to the general result that Maple returned. Although we can’t calculate it directly, We can use limits to see what its value would be at a=1/3, b=0

>simplify(limit(f(1/3, x), x = 0));

(2/9)*Pi*3^(2/3)/((-(1/3)*I)^(1/3)*GAMMA(2/3)*((1/3)*I)^(1/3))

>evalf(%)

evalf(%);
2.230707053+0.1348486379e-9*I

The symbolic expression looks different and we’ve picked up an imaginary part that’s possibly numerical noise. Let’s use more numerical precision to see if we can make the imaginary part go away. 100 digits should do it

>evalf[100](%%);

2.230707051824495741427486519543450239771293806030489125938415383976032081571780278667202004477941904
-0.855678513686459467847075286617333072231119978694352241387724335424279116026690601018858453303153701e-100*I

Well, more precision has made the imaginary part smaller but it’s still there. No matter how much precision I use, it’s always there…getting smaller and smaller as I ramp up the level of precision.

What’s going on?

All I’m doing here is playing around with the same problem in two packages.  Does anyone have any further insights into some of the issues raised?

I am extremely short sighted and have a contact lens prescription of the order of -10 dioptres or so.  Recently, I was trying to explain to my wife exactly how I see the world when I am not wearing my contact lenses or glasses and the best I could do was ‘VERY blurry!’

This got me thinking.  Would it be possible to code something up that took an input picture and a distance from my eyes and return an image that would show how it looks to my unaided eyes?

Would anyone like to give this a go?  How hard could it be?