When you say "columns of different empirical distributions," are you implying that you are working with nominal and/or ordinal as opposed to continuous data?
If so, I wonder if you could recode each column of each distribution into its own variable and then use fuzzy c-means clustering (check out either the -fanny- or -cmeans- package, if you are an R user). In this case, you would compute a distance matrix--presumably with either a Euclidean, squared Euclidean, or Manhattan metric--where the rows are your observations and your columns are you new variables (where, again, each column corresponds to a column from your empirical distributions). This would actually cluster your observations, but I wonder if you could then identify the variables (columns) that clustered observations tend to be distributed over in similar fashions.
Another option would be something like multiple correspondence analysis (package -ca- in R is a good one). This is really more of a factor analysis than it is a cluster analysis insofar as the goal is to reduce the relationships between row and column elements and project them into a shared low-dimensional space rather than to place the elements into any sort of agglomerative dendogram. However, with correspondence analysis, the focus is on variable relationships rather than observation relationships. In this case, you would produce a matrix where your rows and columns are your variables (the recoded "columns" from your original data mentioned before). Run the mjca function from the -ca- package on this matrix to compute the necessary eigenvalue decomposition, and you can see how your column elements "hang together" based on their similar distributions through the row elements and vice versa.
