Using BLAS and LAPACK on Linux

I've written two other related posts that you might also want to read:

Using an NVIDIA GPU with CUDA For Numerical Linear Algebra in Linux

Trying out CSDP on a Xeon Phi Processor

Introduction

Many software packages for scientific computing spend most of their time performing numerical linear algebra operations including matrix-vector and matrix-matrix multiplication, solving linear systems of equations, eigenvalue computations, and matrix factorizations such as the QR and SVD factorizations. These packages often make use of BLAS and LAPACK subroutines to perform these linear algebra operations.

The BLAS library performs basic linear algebra operations such as dot products, vector norms, matrix-vector multiplications, and matrix-matrix operations.  The LAPACK uses functions from the BLAS library to perform higher level computations such as solving linear systems of equations, finding eigenvalues of a matrix, or finding the singular value decomposition of a matrix.  You may also see references to two earlier libraries, LINPACK (for systems of equations) and EISPACK (for eigenvalue problems) that also use BLAS.  They have been effectively superseded by LAPACK.  In most cases, performance optimization at the BLAS level is sufficient.  Some of the libraries mentioned below also include optimized versions of some LAPACK routines, but in general, it's usually sufficient to use an optimized BLAS with the standard LAPACK implementation.

LAPACK and BLAS were originally written in Fortran, and most software that calls BLAS/LAPACK routines uses the original Fortran interface.  For programs written in other languages such as C, this can create complications in interfacing to the libraries.  See the appendix below on calling Fortran BLAS/LAPACK routines from C.  There are alternative interfaces that wrap around the BLAS and LAPACK libraries to provide a more natural C language interface.  The LAPACKE interface provides C interfaces to the LAPACK routines.  The CBLAS interface provides C interfaces to the BLAS routines.

The performance of the software can depend greatly on how well the BLAS operations have been optimized.  There are a number of available implementations of the BLAS with different licenses, performance characteristics, features, and ease of installation and use.  In helping users of my own package, CSDP, I often find that they're struggling to install BLAS and LAPACK to work with CSDP or that they're getting poor performance out of the BLAS that they've used.  

This post discusses the available implementations of the BLAS on Linux and gives some recommendations on how to select a BLAS.  I've written the blog posting from the point of view of someone who is installing software that makes use of the BLAS.  A few appendices at the end of the post discuss technical issues that may be of interest if you're actually writing Fortran or C code that calls BLAS or LAPACK routines.

Factors Affecting the Performance of BLAS

In this section, I'll discuss some important aspects of computer architecture that affect the performance of the BLAS and LAPACK libraries.  Intel and AMD processors based on the x86-64 (aka Intel 64) architecture are most commonly used on desktop and laptop Linux machines, so I'll limit my discussion to these processors. 

It is import to understand that on modern processors it is relatively time-consuming to access a floating point number in main memory- it can take hundreds of nanoseconds to do this, while floating point additions, subtractions, multiplications, and divisions can all be done in a nanosecond or less.
To deal with the slow memory access, modern processors have cache memory.  The cache memory is much smaller than main memory but also much faster.  Processors now have three or even four levels of cache between the processor and main memory.  

The key to obtaining higher performance with cache memory is to make multiple uses of data that has been brought into the cache from main memory before this data is flushed out of the cache.  Algorithms that make a single pass over the data with one or two floating point operations for each number brought from main memory are effectively limited by the speed of main memory and the cache doesn't help.  

The linear algebra operations performed by the BLAS are conveniently categorized into three levels.  The level one operations act on one or two vectors of length n and perform a total of O(n) operations.  Examples include vector dot products, norms, and the saxpy operation in which a multiple of one vector is added to a second vector.  The level two operations involve matrices of and vectors of size n and perform a total of O(n^2) operations.  The most common of these operations is matrix-vector multiplication.  The level three operations involve matrices of size n but perform O(n^3) operations.  The most common level three operation is matrix-matrix multiplication.  

It should be clear that it is impossible to make effective use of cache with level one and level two operations since there are only a small number (one, two, or three) of floating point operations for each number brought in from main memory.  However, the level three operations perform O(n^3) operations on O(n^2) data, so a factor of n reuse of data is possible.  For example, blocked matrix-matrix multiplication algorithms take two large matrices of size n by n and break the matrices into block matrices whose individual blocks are small enough to fit in cache.  The products of blocks are then performed in cache added together and written out to main memory.   Cache aware blocked level 3 operations can be vastly faster than naively coded loops.  LAPACK was designed to make as much use as possible of level 3 operations to improve performance on computers with cache memory.  

An important aspect of modern random access memory is that the memory operates in burst mode.  Rather than reading a single word from memory into the cache, the memory system reads in a "cache line" that might contain 64 or 128 consecutive bytes of memory.   This is very convenient if the algorithm accesses data consecutively but very inefficient when the algorithm jumps around in memory.  Matrices can be stored in memory in "column-major" order with the elements of a column in consecutive positions in memory or "row-major" order with the elements of a row in consecutive positions.   Algorithms used by BLAS and LAPACK are designed to take advantage of the column-major storage used in the Fortran programming language.  See the appendix below for a more detailed discussion of how Fortran and C store matrices.  

In recent years microprocessors have become faster by performing operations in parallel.  There are two main levels at which parallelism has been used to improve performance.  At the lower level are instructions that can perform multiple arithmetic operations simultaneously.  These Single Instruction Multiple Data (SIMD) instructions make it possible to perform 2, 4 or even 8 floating point operations per cycle rather than a single operation.  In successive generations, SSE, SSE2, and AVX instructions have been added as new features to the x86-64 architecture.  Compilers are generally capable of automatically producing code that performs operations in parallel using SIMD instructions.  However, if a code is to be compiled so that it can be distributed in binary form and run on all x86-64 processors then the more recent SIMD instructions cannot be used.  

At a higher level, modern processor chips actually consist of multiple processor cores.  As individual transistors have gotten smaller on integrated circuits it has become possible to cram more and more processor cores onto a single chip.  Some of the latest generation processors have as many as 18 separate processor cores on a single chip.  Some programs are inherently sequential and can only run on one processor core.  In other cases, it is possible to break up the work so that it can be performed in parallel on multiple cores.  Fortunately for us, it is relatively easy to write multithreaded versions of the BLAS functions that can make use of multiple processor cores. 

In order to get high performance, an implementation of the BLAS must make effective use of cache memory, SIMD instructions, and multiple cores.  The software can either be written to work with the most basic x86-64 architecture (so that it will run on all AMD and Intel processors) or written to dynamically select routines that will work well on the current processor, or it can be tuned at compile time for a specific processor.   All three of these approaches are used by some of the BLAS libraries listed below.

Available BLAS and LAPACK Libraries

The following table summarizes the features of some of the most commonly used BLAS and LAPACK libraries.    OpenMP and ILP64 features are discussed in the appendices.

Library	License	BLAS	LAPACK	CBLAS	LAPACKE	Threaded	OpenMP Aware	ILP64
NETLIB LAPACK	BSD		X		X
NETLIB BLAS	BSD	X		X
ATLAS	BSD	X	Partial	X		X
OpenBLAS	BSD	X	X	X	X	X	X	X
MKL	Prop	X	X	X	X	X	X	X
ACML	Prop	X	X			X	X	X

NETLIB LAPACK

The LAPACK project is hosted at the NETLIB web site along with many other packages for numerical computing.  The current version of the library is 3.7.0.   The LAPACK library is written in Fortran, but a C language interface, LAPACKE, is also distributed with LAPACK.  See the LAPACK package in Ubuntu 16.04.

NETLIB Reference BLAS

The reference implementation of the BLAS is a single threaded Fortran library that makes no use of parallel processing and doesn't take advantage of cache memory.  Some performance improvement is possible by using compiler optimizations and SIMD instructions.  The Ubuntu package libblas3 is a generic x86-64 version of the reference BLAS.  It's easy to install but performs poorly in practice.

ATLAS

ATLAS stands for Automatically Tuned Linear Algebra Subroutines.  The ATLAS library comes with a very complicated build procedure that runs timing tests to tune the software to select block size parameters and algorithms that work best on the current machine.  The build process is time-consuming and requires that CPU throttling (commonly used for power management and cooling) be turned off.   Once you've worked through this complicated build process, the library general performs well.  ATLAS includes a few of the LAPACK routines, but for anything complicated, you'll need the regular LAPACK library.

Ubuntu has a generic x86-64 version of ATLAS, but its performance is typically poor in comparison with a version built from source because the generic binary package hasn't been properly tuned.

OpenBLAS

The OpenBLAS library can be compiled from source and automatically configures the block sizes and algorithms for your machine if the microarchitecture is on its supported list.  For newer processors, it's generally OK to use the most recent supported architecture.  For example, I built my copy of OpenBLAS on a Kaby Lake processor using the configuration for a Haswell processor and got good performance.  OpenBLAS includes the cblas interface as well as LAPACK and LAPACKE routines. OpenBLAS works particularly well with main programs that are compiled with OpenMP parallelism.

Ubuntu has a generic x86-64 version of OpenBLAS, but much like ATLAS and the reference BLAS, this package performs poorly because it hasn't been tuned.  (I'm told that this may have been fixed recently with a DYNAMIC_ARCH version of the openblas package.)

OpenBLAS can be built with a DYNAMIC_ARCH flag so that it will dynamically select code optimized for the particular architecture of the machine that it is running on.  A variety of Intel and AMD processors are supported. Executable code built with this version of the library will be quite large, but the resulting code will typically run well on both older and newer systems.  

Intel MKL

The Intel MKL library includes BLAS, LAPACK, and other useful mathematical subroutines.  The cblas and LAPACKE interfaces are also included.  The BLAS is highly optimized and dynamically configures itself on a range of Intel processors.  The library doesn't work well on AMD processors.  Intel has a number of ways of getting free (as in beer) licenses to use MKL.  It is not open source software, however.

AMD ACML

ACML is roughly comparable to Intel's MKL in its functionality but optimized for use on AMD processors.  The library doesn't work well on Intel processors.  The ACML software is proprietary but can be downloaded for free (as in beer.)   AMD has announced the end of support for the proprietary ACML and is moving to support open source high-performance computing projects including the BLIS effort to replace BLAS.  (BLIS is pretty new and complicated to build, so I haven't tried it yet.)

Recommendations

If you simply need to run software that requires the BLAS and aren't concerned about performance, then I would recommend installing a precompiled OpenBLAS package or perhaps a precompiled version of the reference BLAS.  If performance is of only moderate concern, then a precompiled version of OpenBLAS may be suitable.  I would generally not pick ATLAS at this point because of the difficulty of installing it.  If performance is important then I would compile OpenBLAS from source.  The Intel MKL and AMD ACML libraries provide performance similar to (and perhaps better than) OpenBLAS, but I wouldn't bother with MKL or ACML unless you also need access to other features of these libraries beyond BLAS.

For example, I recently installed my CSDP software on a newly purchased laptop with a Kaby Lake (latest generation Intel) processor.  The Ubuntu OpenBLAS binary package performed poorly.  I found that a version of OpenBLAS compiled from source was at least four times faster than the Ubuntu OpenBLAS binary package, the Ubuntu ATLAS package, and the Ubuntu reference BLAS package.  ACML performed poorly on this Intel machine.  MKL was roughly similar in performance to the version of OpenBLAS that I built from source.

Appendix: Fortran and C array storage

In Fortran, an m by n matrix A is stored in a data structure called an array.  To give flexibility in adding rows to a matrix, Fortran and the BLAS/LAPACK interfaces allow for an array A to have LDA > m rows.  (LDA is an acronym for “leading dimension of A.”)

It is necessary to map the two-dimensional array of entries of A onto the one-dimensional linear layout of memory.  In the Fortran scheme, the array is stored in “column major” order. That is, the entries in the first column are followed by the second column, and so on.  In general, the A(i,j) element of the array A appears in position i+(j-1)*LDA in memory.  

For example, a 3x3 matrix might be stored in a 5x3 array as

1	A(1,1)
2	A(2,1)
3	A(3,1)
4	X
5	X
6	A(1,2)
7	A(2,2)
8	A(3,2)
9	X
10	X
11	A(1,3)
12	A(2,3)
13	A(3,3)
14	X
15	X

In calling Fortran routines from the BLAS and LAPACK libraries, the programmer must pass the array A together with the matrix size parameters m and n and the leading dimension LDA.

C actually has two different schemes for storing arrays.  Dynamically allocated arrays consist of a vector of pointers to 1-D arrays which represent the individual rows of the matrix.  This style of arrays is not supported by the cblas and LAPACKE interfaces.  Statically allocated two-dimensional arrays are stored in a scheme very similar to the Fortran arrays, except that the array is stored in “row major” order, with all of the entries in row one proceeding row two, and so on. Another less important difference is that C array subscripts start at 0 rather than 1.

A C programmer has two options for interfacing to BLAS and LAPACK.  The older approach is to implement Fortran array indexing within 1-D C arrays and then call the standard Fortran BLAS/LAPACK routines.  This approach is used by many software packages. The alternative approach is to use LAPACKE and cblas if they are available. The LAPACKE and cblas interfaces work with matrices that are stored in either row-major or column-major order.  

Appendix: Calling Fortran subroutines from C

Calling Fortran subroutines from C is not completely standardized.  Different Fortran compilers use different naming and parameter passing conventions, so it is sometimes necessary to adjust C code to work with a different Fortran compiler.  

In Fortran, subroutine names are case-insensitive (and often converted into lower-case by the compiler), while in C, they are case sensitive.  Many Fortran compilers insert an underscore ("_") at the end of subroutine name.  

In Fortran, parameters are passed by reference (the address of the parameter) rather than value, while in C parameters are typically passed by value.  Thus when you call a Fortran subroutine from C it is necessary to pass a pointer to each of the parameters.    

In C, variable length strengths have a null byte (0) to indicate the end of the string.  Some Fortran compilers include additional hidden string length parameters to specify the length of each of the strings that are being passed to the subroutine.  

Consider the Fortran call to an LAPACK routine

   CALL DPOTRF('U',N,A,LDA,INFO)

where the string 'U' tells LAPACK to put its answer in the upper triangular part of A, N is the dimension of the matrix A, LDA is the leading dimension of the array A (not necessarily equal to N!), and INFO returns a result code.

In C, the most common way to write this would be

 dpotrf_("U",&n,A,&lda,&info);

However, depending on the Fortran compiler, it could be

  dpotrf("U",&n,A,&lda,&info);

or

  DPOTRF("U",&n,A,&lda,&info);

or

  DPOTRF_("U",&n,A,&lda,&info);

or

  dpotrf_("U",&n,A,&lda,&info,1);

In this last example, the parameter "1" specifies the length of the string "U".


Appendix: Pthreads and OpenMP threading

Many of the BLAS libraries are "threaded", meaning that the library will automatically spawn several threads of execution on different processor cores to compute (e.g.) a matrix-matrix multiplication.  A program that uses BLAS/LAPACK libraries might itself run single threaded, or it might itself be multithreaded. Typically, it's best for a program to use as many threads as there are available processor cores.  Running more threads than cores can hurt performance as the computer must spend time switching back and forth between threads.  A common performance problem occurs when a multithreaded main program has several threads that each call threaded BLAS routines.  E.g. on a four core machine, a program might use 4 threads, but these 4 threads can each call a BLAS routine that starts 4 threads, and suddenly there are 16 active threads.

There are two very commonly used threaded programming models.  The pthreads (posix threads) model and OpenMP threading.  Intel also has its "threaded building blocks", although these are not very widely used in open source software.   There is no standard way in pthreads for a process to indicate that subroutines shouldn't create new threads and thus no effective way to avoid the performance problem mentioned in the previous paragraph.  However, OpenMP does have a standard mechanism for doing this.  Some BLAS libraries are "OpenMP aware" and will not create additional threads if the calling program doesn't want additional threads.  This is a very useful feature.

Appendix: I32LP64 vs. ILP64

On 64 bit Linux systems, C uses 32 bits to represent integers, 64 bits to represent long integers, and 64 bits for pointers.  Similarly, Fortran integers are 32 bits.  This programming model is referred to as I32LP64.  An alternative model, ILP64, uses 64 bits to represent integers.  There is a potential issue with using the I32LP64 model in BLAS/CBLAS and LAPACK/LAPACKE, since it is possible in some circumstances that a vector might be so long that we can't represent its length with a 32-bit integer.  The largest number that can be represented by a signed 32-bit integer is about 2 billion, but BLAS and LAPACK are often used on machines with hundreds of gigabytes of RAM, so this actually can occur in some cases.  Some implementations of the BLAS and LAPACK libraries have versions that use or can be compiled to use long int 64-bit integer parameters.

Appendix: CPU Frequency Scaling

Modern processors from Intel and AMD have features to adjust the processor clock rate as needed to save power and prevent the processor from overheating.  The processor cores can also go into a faster than normal (but not sustainable for long periods) "turbo" mode.   For example, the cores on my desktop computer can run at clock rates from 1.0 Ghz (not loaded) to 3.7 Ghz (the nominal clock rate of the processor) and up to 4.0 Ghz (in turbo mode.)  Rather than operating at a fixed, high, clock rate, processor cores are usually kept at a low clock rate.  When the system has to start performing computations, the clock rate is adjusted upward by the operating system.   However, if the processor cores start to overheat, the system will turn the clock rate down.   

As a system administrator, you can control the minimum and maximum clock speeds and the "governor" procedure that determines how the system responds to increased mode.  The default governor, "ondemand" increases the clock rate in response to high system loads, but it can take some time to get up to full speed.  The "powersave" mode is useful for laptops running on battery power but doesn't give very high performance.  There is also a "performance" mode which is the mode should typically be used for scientific computing.  The cpufreq package has command line tools that can be used to configure CPU frequency scaling.    I keep my systems in "performance" mode.  I've configured my laptops for "powersave" while unplugged and "performance" while plugged in and charging. 

Appendix: Hyperthreading (Simultaneous Multithreading)

Some AMD and Intel processors are capable of "simultaneous multithreading (SMT)" (Intel marketing uses the term "hyperthreading" for this.)  This means that a single core can simultaneously execute two (or more) threads.  The processor core has additional registers so that it can quickly switch from one thread to another.  For example, when one thread of execution is stuck waiting for data from memory, the other thread can take over immediately without any of the overhead associated with switching tasks at the operating system level.  For number crunching, hyperthreading typically provides relatively little benefit, and in some cases, codes even run slower with hyperthreading turned on.  It can also be confusing that a 4-core processor with 2-way hyperthreading appears to Linux as an 8-core processor.  You should experiment to see how much if any performance improvement you get by running more threads than processor cores.  

Appendix: The TLB and Transparent Hugepages

Modern processors make use of "virtual memory", which allows the operating system to move pages of code and data in and out of main memory and store them on disk.  This makes it possible to have multiple programs running and using more "virtual memory" than the real memory in the system.   Because pages can end up being moved to different real memory locations as they are paged in and out of main memory, it is necessary to translate virtual addresses into real physical addresses.  The operating system maintains an address translation table that maps each virtual memory page's virtual addresses into a corresponding physical address (or indicates that a page is not in main memory and must be brought in from disk before being used.)

A special cache called the translation lookaside buffer (TLB) stores the address translation for a small number (at most a couple of thousand) of the most recently accessed pages.  If a page's information isn't in the TLB, then the operating system executes code to bring that entry from main memory into the TLB.   The TLB handles "hits" instantly,  while TLB "misses" cost a considerable amount of time (a microsecond or more.)

The Linux kernel was originally designed to use only 4K byte pages.  More recently, support for hugepages (of 2 megabytes and also 1 gigabyte for the latest processors) was added to the kernel.  Programs can explicitly allocate the larger pages, or a Linux feature called "transparent hugepages" can be used to dynamically identify situations where small 4K byte pages can be combined into larger hugepages.

For numerical linear algebra codes involving multi-gigabyte matrices, the older 4K byte pages are very problematic.  A 4-gigabyte array uses one million 4K byte pages, far more than the capacity of the TLB.  This results in many performance-sapping TLB misses.  With two megabyte huge pages, the 4 gigabyte array uses a much more manageable two thousand hugepages.  For some other workloads, particularly databases and web servers, that rapidly allocate and deallocate storage, the transparent hugepages feature actually hurts performance.  The decision to turn on or off transparent hugepages depends on the workload of a particular computer.  Some Linux distributions enable the feature by default while others disable it.

You can check to see whether transparent hugepages are enabled on your Linux machine with

cat /sys/kernel/mm/transparent_hugepage/enabled

The answer will either be "always", "madvise" (if a program has to ask for huge pages) or "never."  If desired, the transparent hugepage feature can be turned on by an administrator by

sudo hugeadm --thb-always

Using transparent hugepages can improve the performance of codes that use BLAS and LAPACK on large arrays.  For example, I typically see about a 10% performance improvement with hugepages on my laptop computer.  

Appendix: The Alternatives system on Debian and Ubuntu Linux

Debian Linux (and by inheritance, Ubuntu) has a feature that makes it possible to have multiple versions of BLAS, LAPACK, and other libraries installed and to automatically adjust symlinks as packages are added to the system that implement these libraries.  For example, the standard location for the version 3.x LAPACK library is /usr/lib/liblapack.so.3.  This is actually a symlink to /etc/alternatives/liblapack.so.3, which itself is a symlink to the currently active version of LAPACK.  

You can query the available alternatives for a package and use the update-alternatives tool to update the configuration as desired.  For example, 

sudo update-alternatives --query liblapack.so.3

will tell you the current configuration for liblapack.so.3.  You can pick from among the installed versions of LAPACK with 

sudo update-alternatives --config liblapack.so.3

If you have your own version of openblas in /opt/openblas/libopenblas.so, you can make it the active version of liblapack.so.3 by 

sudo update-alternatives --install /usr/lib/liblapack.so.3 liblapack.so.3 /opt/openblas/libopenblas.so 99

where the priority "99" should be large enough to keep other versions from taking over as liblapack.so.3.  


Further Reading

For general background on computer architecture including caching, the translation lookaside buffer, and other features that are important to performance, I highly recommend the book:

R. E. Bryant and D. R. O'Hallaron.  Computer Systems: A Programmer's Perspective, 3/E, 2015.   An earlier draft edition is available for free online.

For background on numerical linear algebra, algorithms that work by rows versus columns, and blocked algorithms, I highly recommend:

D. S. Watkins.  Fundamentals of Matrix Computations, 3rd ed.

Scraping Citation Counts from a Google Scholar Profile

Google Scholar has no published API.  This means that if you want to extract information from someone's public Google Scholar profile, you have to write a script to scrape the web page.

I was recently asked to come up with a way to extract citation counts by year from a Google Scholar profile. This information is visible in a chart on the profile page, and it wasn't too hard to write a Python script to extract it. You can see my script at

http://www.nmt.edu/~borchers/scrapescholarcitecounts.py

Note that the citation counts are actually plotted as percentages of the citation counts in the highest year, and these percentages are rounded to one digit after the decimal point.  e.g. if an author's maximum citation year had 70 citations and a particular year had only 3 citations, then this is encoded in the table as 4.3%.  The script rounds to the nearest integer, but this can easily be off when the total number of citations is high.  

Similar plots are produced by Thompson's "MyResearcherID" profile, but in that case the chart is produced by a process that is opaque- you can get at the produced GIF, but in order to extract counts by year you'd have to run this graph through a "graph digitizer" tool.

It's a sad commentary on the state of the world wide web that information like this is presented in web pages in formats that are very hard to decode.  

Review of Coursera's "Computing for Data Analysis" MOOC

I signed up for this course primarily because I wanted to experience a Coursera massively open online course (MOOC), and also because I wanted to improve upon my very rudimentary knowledge of the R programming language.  For those of you too impatient to read the entire review the summary is that I did learn a few things from the course and I'm glad that I took the time to participate, but there were many things about the course that really could be improved.

First, a comment about the word "open" in the acronym MOOC.  Although some of the first MOOCs really were open in the sense that the course content was made available under open source licenses, the Coursera course content is not made available under any open license.  This makes it effectively impossible for students and other instructors to build upon the MOOC and improve it.  Many of my criticisms of Computing for Data Analysis could be very effectively addressed by incremental improvements made in the typical fashion of source software projects.

It is interesting that there are already many open resources for learning R programming in the form of online tutorials and notes. These clearly aren't as interactive and engaging as the Coursera course, but it seems likely that the authors who produced these materials have enough knowledge of the subject so that with the support of instructional designers and producers they could also produce a MOOC on R programming.  Thus the special sauce is clearly not in Roger Peng's expertise in R- lots of instructors have that background.

I also don't believe that the Coursera platform is a unique distinguishing feature- Coursera and its competitors have all been able to produce MOOC platforms that can scale to handling hundreds of thousands of students with very little difficulty.  The technologies involved in these platforms, including lecture videos with subtitles and indexing, online multiple choice tests, banks of test questions, algorithmically generated questions, and even more sophisticated adaptive instruction techniques are already available in many computer based instruction systems.

Many universities already make use of the same textbooks and instructional software- they're not competing with each other by controlling access to the content that they're teaching.  Rather, they're competing with each in other areas, including the quality of interaction between students and faculty, the scholarly reputation of faculty, and various aspects of the college experience unrelated to academics.  In fact, faculty are often encouraged to publish textbooks because this helps to bolster the reputation of an institution.  If universities don't want to compete in the area of content, but they're concerned about Coursera dominating the world of MOOCs, then they have an easy alternative.  They can cooperatively develop truly open content in the form of Open Educational Resources (OER).

OER also offers an escape route for faculty who fear that they may be reduced to nothing more than tutors and graders for courses in which an academic star is featured in recorded lectures.  With OER, instructors have the opportunity to incorporate previous work by other instructors into their own courses, but they also have an opportunity to bring their own work into the mix.  Everyone gets to participate, and instructors can focus on finding the best ways to present the
material that they are most interested in.

For these reasons, I think that OER represent a real and significant alternative to the commercial MOOCs.

Now on to the Computing for Data Analysis (hereafter I'll just use the acronym CFDA)  course itself.

One obvious issue in constructing a MOOC is whether to set the course up in self paced format or to run it on a schedule.  CFDA was run on a four week schedule, with weekly quizzes and programming assignments.  Given the way that the video lectures, quizzes, and programming assignments were presented, there was no technical reason that the course couldn't be presented in a self paced format with students working at their own pace.  The scheduled format did make it easier for participants to discuss the course in the online forums associated with the course.  However, it also caused problems for students who had to travel or were otherwise busy at some point during the four weeks of the course.  To help with this, students were given loose deadlines and a small number of extra days that they could spend on assignments.  I personally found the scheduled format to be overly restrictive, and I would have preferred to see the course presented in a self paced format.

The lectures in this course were presented as online videos of up to about 1/2 hour length.  The videos consisted of prepared slides with the lecturers voice in the soundtrack.  To keep the videos from becoming too boring, there were multiple choice questions interspersed throughout the videos.  There were also some activities for the viewer to try at the R command line.  The production values were comparable to what I regularly see in lecture capture videos from conventional courses-  these were not highly polished videos.

In most conventional courses, lectures are at least loosely tied to a textbook.  Although there are many open source texts on R and also many commercially available textbooks, CFDA was not tied to any particular text.  Using a  commercial text would obviously have been problematic in a "free" course.  Using one of the open source resources would have been equally problematic given Coursera's closed source approach.

Not having a reference is problematic in a course about a programming language, since anyone learning a programming language is going to have to frequently look up the details of syntax and library functions.  It was quite irritating to me in doing the programming assignments that I frequently had a vague memory of some function from the lecture but didn't remember it in enough detail.  I then had the unpleasant choice of going back through the video or going to some other resource to find that information.

I was disappointed by the lack of depth in the lecture material.  This was a course in R programming rather than statistics, so any in-depth discussion of statistical issues would have been inappropriate, but I was surprised by how little discussion there was of the structure and semantics of the R programming language.  For example, R uses a fairly unusual scheme of lexical scoping that can be particularly useful in statistical work.  This was mentioned briefly in one  lecture but not fully explained.

The quizzes were presented in multiple choice format, with the student given  the opportunity to retake a quiz up to 3 times to get a higher score.  Some of the  questions were written in such a way that it was easy to guess the correct answer without actually understanding the question.  On retakes, the same questions reappeared, sometimes with the order of the answers randomized.  The system did not appear to be using algorithmic question generation or even a large library of similar questions.   Although multiple choice questions like these are easy to implement in a MOOC, they simply don't have the depth of questions that require a written answer or actual coding.

The programming assignments were automatically graded- the student would simply run test scripts that tested his code on various example inputs and then upload the resulting output for checking.  Although this kind of automated grading of programming assignments can be quite rigorous if its done with carefully designed tests, the test cases used in this course seemed to be pretty easy.  None of my codes ever lost a point.

In my opinion, this kind of automated grading of programming assignments often leads to poor programming practices, since it pushes students to debug until they "get the right answer", rather than developing a program that know to be correct by analysis.  Furthermore, there's no opportunity to give the student feedback on the proper style of programming in the particular language.  The
programs that I wrote for CFDA were not "pretty", and probably would have come across to an expert R programmer as being non-idiomatic R code.  All that mattered within the course is that they produced the correct answers on the test cases.

Although the online discussion forums were billed as one of the most important features of the course, I was very disappointed by them.  The way that the assignments were set up, students were focused on completing particular tasks rather than on developing a broader understanding of the R system.  As a result, much of the discussion in the forums was at the level of "How do I do in R?"  I found that I was able to more efficiently answer my questions of this sort by using other online resources on R and a bit of Google searching.

Although I think this course could have been much better, I certainly did learn some things about R programming that I wanted to pick up, and overall I'm glad that I took the time to participate in the course.

Thoughts on Tablets

I've recently spent a couple of weeks trying out an iPad 2 followed by a couple of weeks testing out a Motorola Xoom tablet.  Despite my dislike for Apple's technology and business models, I have to say that the iPad 2 is definitely a much more functional device at this point in time.  However, neither device was really satisfactory from my point of view.

If I wanted a device of tablet size for web browsing, email, calendar, video calling, and ebook/PDF reading, then I'd buy an iPad 2- it does all of this just fine.  Now that IOS 5 goes tetherless, I could use an iPad without needing to keep a Windows machine around to run  iTunes- this is a big plus.

The Xoom also did web browsing, email, calendar, etc.  However, it just wasn't as well polished and not nearly as fast- the performance of the Xoom was quite disappointing.  I had hoped that Google integration on the Xoom would be better than on the iPad, but it really wasn't any better.

I've given both of these borrowed devices back to the people that I loaned them from, without any regrets.

I'm quite happy with my current Android phone.  It's always in my pocket, and does almost as much as an iPad or Android tablet while being vastly more convenient to carry around.  Until the tablets offer a lot more functionality, I don't think I'm going to acquire one.

The biggest things I want that I'm not seeing in the current generation tablets:

1. A stylus for handwriting with reasonable resolution (this is particularly important because I hand write a lot of mathematics, particularly in my teaching.)  

2. Flash and Java applets.  Neither are available on the iPad, and Java applets don't run on Android.  Android support for flash is somewhat limited.  Of course we hope that HTML 5 will kill off flash once and for all, but there doesn't seem to be any reasonable replacement for Java applets on the horizon.  Unfortunately, lots of good educational software uses these technologies.

3. Wireless display to computer projectors (for presentations and teaching.)  Apple has a new offering in this direction, but I'm looking for a more open standard, such as Wi-Di.

For this list of features it appears that Windows Tablet PC's (an old class of machines dating back 10 years!) are still the best way to go.  These are heavier devices than the Android and IOS tablets, but they feature much more powerful processors and more interface options.  A few years ago these were two thousand dollar devices, but there are some reasonably good models out now in the 1,200 dollar price range.  I was particularly impressed with how well Microsoft's OneNote software works on a tablet PC.

Convex Optimization 2.0

My sabbatical leave last fall at the Institute for Pure and Applied Mathematics (IPAM) at UCLA was a wonderful opportunity to get up to speed on some recent developments in convex optimization and to reassess where the field is going. I recently gave a talk to our math department seminar with the somewhat provacative title of Convex Optimization 2.0.

This posting is a brief summary of what I talked about, along with some links to useful resources related to the talk.

I'll start with a quote from R. Tyrrell Rockafellar in 1993.

In fact, the great watershed in optimization isn't between linearity and nonlinearity, but convexity and nonconvexity.

R. Tyrrell Rockafellar (1993)

At the time, this was a somewhat surprising statement. Through the 60's, 70's, and into the 80's, most people working in the field thought of mathematical optimization being divided into two parts- linear programming and nonlinear programming. The central mathematical ideas of linear programming were combinatorial in nature, while the mathematics of nonlinear programming was all about analysis of smooth functions and applications of Newton's method.

What's amazing to me about this quote is that today, most of the people that I know who are working in optimization would agree with Rockafellar. The point of my talk and this blog post is to tell the story of how we got to this point.

George Dantzig's simplex method was the preferred method for solving linear programming problems from the 1950's through the mid 1980's. An enormous amount of research and software development effort (and lots of CPU cycles) were invested in the simplex method.

Methods for nonlinear programming are often divided into first order methods that use only first derivative information and second order methods that use second derivative information to obtain faster convergence. Of course second order methods also require storage that grows quadratically with the size of the problem. Thus first order methods are often preferred for very large scale problems, while second order methods are more commonly used for small to medium sized problems. There are also methods such as the Limited Memory BFGS Method that try to fit into the space between first and second order methods to get faster convergence than the first order methods without the memory requirements of second order methods.

When primal-dual interior point methods for linear programming appeared on the scene in the mid 1980's, it was surprising that ideas from nonlinear programming could actually be useful in linear programming. In less than a decade, these new interior point methods evolved to the point that they were often competitive with the classical simplex method, particularly on very large and degenerate linear programming problems. Interior point methods were also being applied to nonlinear programming problems, further blurring the distinction between linear and nonlinear programming.

Convex optimization problems are problems in which the feasible region is a convex set and the objective function to be minimized is a convex function. Some researchers, including Rockafellar, had worked on methods for solving nonlinear optimization problems that were convex but not necessarily smooth. The huge advantage of this class of problems is that any local minimum of a convex optimization problem is also a global minimum. On the other hand, without smoothness it is impossible to directly apply methods like Newton's method that depend on the differentiability of the objective and constraint functions. Methods for nonsmooth convex optimization (or "nondifferentiable optimization" as it was sometimes called) were quite limited.

After the development of interior point methods for linear programming, it was quickly discovered that the primal-dual interior point method for LP could also be extended to handle other kinds of convex optimization problems. For example, interior point methods can be used to solve convex quadratically constrained programming problems (CQP).

Of particular interest were second order cone programming problems which involve constraints of the form

\(
\|Ax-b \|_{2} \leq c^{T}x+d
\)

and semidefinite programming problems in which the variable $X$ is a symmetric matrix that is constrained to be positive semidefinite.

Looking at this situation from the point of view of our ability to solve various kinds of convex optimization problems, we now had a heirarchy

\(
LP \subset CQP \subset SOCP \subset SDP \subset CP.
\)

with primal-dual interior point methods able to efficiently and reliably solve problems up to the level of SDP.

Although SOCP and SDP might seem to be extremely weird problem classes, it turned out that in practice, many convex optimization problems could be reformulated into SOCP or SDP form. This led to a paradigm that I'll call Convex Optimization 1.0. Start with a convex optimization problem, reformulate as an LP, CQP, SOCP, or SDP, and then use a primal--dual interior point method to solve the problem. This approach is beautifully laid out in the 2003 textbook, Convex Optimzation, by Stephen Boyd and Lieven Vandenberghe.

In fact, the reformulation techniques used to put a problem into standard form are so simple that they can actually be built into software. For example, Michael Grant and Stephen Boyd developed a MATLAB package called CVX that automatically reformulates convex optimization problems as LP's, SOCP's, or SDP's, and then solves them using a general purpose LP/SOCP/SDP solver routine. Some other important software packages that also implement this approach include YALMIP, ROME for robust optimization, and GLOPTIPOLY and SOSTOOLS for polynomial optimization problems.

As interest in convex optimization grew, one important problem with primal-dual interior point methods became apparent. The primal-dual interior point method is a second order method, with storage requirements that grow as the square of the problem size. Users of interior point codes were often disappointed that the problems they wanted to solve simply couldn't be solved because of this storage requirement. The need for more storage efficient methods was particulary important in the burgeoning field of compressive sensing, where researchers wanted to solve problems involving millions of variables.

This led to a revival of interest in first order methods for convex optimization problems. Although these methods had originally been developed in the 1970's, this was an area of research that had been basically dead for a long time.

The classic algorithm for smooth convex minimization problems is the projected gradient descent method. Given a convex function $f(x)$, and a convex feasible set $C$, and a function for computing the projection of a point onto the nearest point in $C$, the projected gradient descent method simply alternates between moving downhill and projecting back onto the set $C$. If the problem is unconstrained, the $P_C$ is trivial. The projected gradient descent method follows the iteration:

\(
x^{k+1}=P_{C}(x^{k}-\alpha_{k}\nabla f(x^{k}))
\)

The step lengths $\alpha_{k}$ can be chosen in many different ways to insure convergence. One simple rule is $\alpha_{k}=1/k$.

For nonsmooth convex optimization problems, the projected gradient descent method can be extended by using subgradients instead of the gradient vector. A vector $g^{k}$ is a subgradient of $f$ at the point $x_{k}$, if

\(
f(x) \geq f(x^{k})+g^{T}(x-x^{k})\;\; \mbox{for all} \; x.
\)

The convergence rate of these methods is sublinear. However, it can be useful to analyze how many iterations are required to obtain a solution $x^{k}$ with $f(x^{k})-f(x^{*})<\epsilon$, where $x^{*}$ is an optimal solution. The classical (projected) (sub)gradient descent method takes $O(1/\epsilon)$ iterations on smooth problems and $O(1/\epsilon^{2})$ iterations on nonsmooth problems.

Nesterov (1983) gave an optimal first order method for smooth convex optimization problems with Lipschitz continuous gradients. This accelerated first order method uses previous gradients in a recursive formula for computing the next step. This effectively provides the algorithm with some information about the curvature of the function. This algorithm requires $O(\sqrt{L/\epsilon})$ iterations.

In 2005, Nesterov returned to this earlier algorithm and showed how a nonsmooth convex optimization problem could be approximated by a smooth convex optimization problem with Lipschitz continuous gradient. Nesterov's smoothed problem had $L=O(1/\sqrt{\epsilon})$, so the accelerated version of gradient descent on the smoothed problem took $O(1/\epsilon)$ iterations. Since 2005 there has been an explosion of interest in methods for various convex optimization problems that use this idea of accelerated gradient descent on a smoothed version of a nonsmooth convex optimization problem.

This approach has been particularly successful in the field of compressive sensing. Here the goal is to find a sparse (few nonzero entries) solution to an underdetermined linear system of equations $Ax=b$. Minimizing the number of nonzero entries in the solution is actually an NP-Hard combinatorial optimization problem. However, it turns out (and this is another story worthy of its own blog posting) that this can be done very effectively in practice by solving a basis pursuit denoising problem .

\(
\min \| x \|_{1}
\)
subject to
\(
\| Ax - b \|_{2} \leq \sigma .
\)

Although these problems can be solved in the Convex Optimization 1.0 paradigm by a general purpose SDP solver, in practice, specialized first order methods are vastly more efficient. Stephen Becker has produced a very useful collection of resources on software for BPDN and related problems.

Now we have the beginings of a new paradigm that I'll refer to as Convex Optimization 2.0. Use an accelerated first order method with smoothing to solve your convex optimization problem.

So far, most of the work on first order methods has been specialized to particular problems. An interesting recent development is the release of Templates for First Order Conic Solvers (TFOCS) by Stephen Becker, Emmanuel J. Candès and Michael Grant. TFOCS deals with these problems by formulating them in conic form (using for example SOCP constraints), dualizing, and then applying a first order method to the dual problem.

Although first order methods have worked out very well in many applications, we haven't yet seen these methods develop to the point that a general purpose first order method for LP, SOCP, and SDP can take the place of a primal-dual interior point method for general problems. The semidefinite programming case seems to be particularly hard. However, I'm hopeful that we may soon see a resurgence in interest in convex optimization as first order methods become more broadly applicable.

LaTeX in blogger postings

I've spent some time this evening setting this blog up so that I can include mathematical formulas like

\(
y=\sin(x)
\)

This is done using the wonderful MathJax package, which allows me to type LaTeX formulas into a blog posting and have them automatically processed by a javascript program that converts the LaTeX for display.

This should work with Internet Explorer, Firefox, Safari, and Google Chrome. I've only tested it with Firefox and Chrome. If you see this and and can't read the above formula, please let me know what browser you're using.

If you're curious about how to do this blogger, please send me an email- the required HTML code is a bit of a pain to describe in a blog posting (where it would end up getting interpreted!)