New implementation of cov/cor with extended API

[lindahua]

This PR is made based on #5249.

• Here is the extended API:

cov(x[, corrected=true, vardim=1, zeromean=false])
cov(x, y[, corrected=true, vardim=1, zeromean=false])
cor(x[, vardim=1, zeromean=false])
cor(x, y[, vardim=1, zeromean=false])

covm(x, xmean[, corrected=true, vardim=1])
covm(x, xmean, y, ymean[, corrected=true, vardim=1])
corm(x, xmean[, vardim=1])
corm(x, xmean, y, ymean[, vardim=1])

Here, the arguments x and y can be either vectors or matrices.

Below is the explanation of keyword arguments:
- corrected: scale the matrix by 1 / (n - 1) if set to true (by default), or 1 / n otherwise. Note that this option does not apply to cor and corm, as it does not affect the resultant correlation values.
- vardim: consider variables in columns and observations in rows if set to 1 (by default), or the other way round, if set to 2.
- zeromean: x and y are considered to have zero means, if set to true. Otherwise, the means will be computed from x and y and subtracted therefrom, if set to false (default).

• The computation is internally delegated to covzm and corzm, which take centered data. These functions further delegate to unscaled_covzm, which is the core implementation. All these internal functions are not exported.
• The functions are thoroughly tested -- I expand the test cases to ensure correctness under all combinations of keyword arguments.

[johnmyleswhite]

This seems like a much better API.

[lindahua]

The PR is ready for review.

[jiahao]

👍

[nalimilan]

Nice!

[StefanKarpinski]

This is a much better API but I really would love to figure out a way to get rid of the *m variants. It was awkward with varm but now there's a whole parallel family, which is worse. Every time there are parallel families of functions like this, it's a major API design smell.

[StefanKarpinski]

You used a new type for weighted vectors in StatsBase – I wonder if something like that is warranted here? I.e. a vector of values with known mean? Is that an approach we can generalize? It's kind of elegant the way it works with multiple dispatch. In that case, it may make sense to move the known mean cases out of Base, but I think that's not unreasonable.

[lindahua]

In current implementation, cor relies on corm, as something like

cov(x) = covm(x, mean(x))

So it makes sense to have covm to live in Base.

Another way to design the API without the need of *m functions may be the following:

cov(x)
cov(x, y)
cov((x, xmean))
cov((x, xmean), (y, ymean))

When people want to additionally provide a known mean, he may use a tuple (x, xmean). This way also uses multiple dispatch, without introducing any new types.

[StefanKarpinski]

I like that tuple approach but I wouldn't mind getting an opinion from @JeffBezanson. I feel like he might have a reason to object that isn't occurring to me.

[lindahua]

If no objection, I will implement the tuple-based interface tomorrow and update the PR.

[lindahua]

Consider this more, with the zeromean keyword argument, it doesn't seem that we need to use covm or other APIs that allow people to input a known mean.

Suppose you already know the mean xmean, you can always write:

cov(x .- xmean; zeromean=true)

or

z = x .- xmean
cov(z; zeromean=true)
# z can be reused later

I am considering just deprecating the *m functions altogether. Opinions?

[nalimilan]

How smart! The only gotcha is that subtracting the mean implies traversing the vector once and allocating space before calling cov, while the *m variants could do the subtraction on the fly if they wanted (though I see that currently you don't do this).

Personally I really don't care about this very subtle potential difference in performance, and the variants could easily be added later if needed.

[lindahua]

Currently, the implementation of covm actually just does the subtraction for the user. The varm does everything in one pass. But the performance difference doesn't seem to be very big.

My proposal: modify this PR a little bit, remove covm and corm, which simply just does x .- xmean for the user, nothing more. The varm and stdm may continue to stay if we don't want to get rid of them for now.

[StefanKarpinski]

How about making the keyword mean instead:

julia> var(x; mean=Base.mean(x)) = (x,mean)
var (generic function with 1 method)

julia> var([1,2,3])
([1,2,3],2.0)

To specify the mean just do something like var(x, mean=0).

[nalimilan]

@StefanKarpinski Isn't it what you proposed here? #5249 (comment) I like this syntax, but it's not clear whether the performance hit due to the keyword argument can be eliminated in the future (cf. the discussion following your comment)..

[lindahua]

I am fine with the mean keyword argument. In particularly, we can do

cov(x; mean=m) = covm(x, m)

Even if the type of mean is not specialized, the only overhead is to resolve the proper covm method in runtime. As the computation in cov is usually heavy, such overhead is generally negligible. In this way, we can have covm as an internal implementation function, which need not be exported.

The problem remains API design, what about cov(x, y). Shall we use xmean and ymean, as?

cov(x; mean=m)
cov(x, y; xmean=mx, ymean=my)

[StefanKarpinski]

How about mean is a triple of the same size as the number of variables?

[nalimilan]

Yeah, the 2-uple sounds like the best solution.

[jiahao]

How about cov(x, y; mean=(mx,my))?

[lindahua]

@jiahao, that's what @StefanKarpinski suggested. I am implementing this. Will be ready soon.

[lindahua]

This is the new API:

cov(x[; vardim=1, corrected=true])    # default: using Base.mean to compute mean
cov(x; mean=0[; vardim=1, corrected=true])  # zero mean
cov(x; mean=m[; vardim=1, corrected=true])   # user provided mean

cov(x, y[; vardim=1, corrected=true])    # default: using Base.mean to compute mean
cov(x, y; mean=0[; vardim=1, corrected=true])  # zero mean
cov(x, y; mean=(mx, my)[; vardim=1, corrected=true])   # user provided means

cor(x[; vardim=1])    # default: using Base.mean to compute mean
cor(x; mean=0[; vardim=1])  # zero mean
cor(x; mean=m[; vardim=1])   # user provided mean

cor(x, y[; vardim=1])    # default: using Base.mean to compute mean
cor(x, y; mean=0[; vardim=1])  # zero mean
cor(x, y; mean=(mx, my)[; vardim=1])   # user provided means

The mean, corrected and vardim keyword arguments are also added to var and std.

[lindahua]

I think it is ready to merge. If no further API change is needed, I will merge this soon.

[nalimilan]

This API looks perfect to me. ;-)

[StefanKarpinski]

My one last (I swear) possible bikeshed is the name of the vardim keyword. Would this be better to call variable? Is vardim clearer or just a weirder name? Sorry for the additional bikeshed. Once we settle on this it would be great to never have to change it ever again :-)

[lindahua]

vardim indicates the dimension of variables. So it has the meaning of dimension.

When vardim = 1, it indicates that the variables are in columns, otherwise (vardim = 2) the variables are in rows.

The name variable does not convey the meaning of dimension. At first glance, variable = 1 seems to mean the variable value is 1.

[StefanKarpinski]

That seems like a reasonable argument. Merge at will!

[toivoh]

Just to continue thd bikeshed :) In numpy, this kind of argument is
invariably called axis. I think that we should be similarly consistent.

[StefanKarpinski]

I dunno. I don't find the name axis to be enlightening at all.

[toivoh]

We're talking about which dimension to reduce along. Although I like
axis, I primarily think that we should be consistent with reduction
functions and others that take the same kind of argument. I think that
dim would be fine; I don't see that vardim adds anything beyond that.

[StefanKarpinski]

Oh, I see what you're saying – not that we should use axis but that we should be consistent.

[lindahua]

There are two kinds of dimensions. The dimension of each sample/observation, or the dimension of each variable. I am just being explicit here.

[StefanKarpinski]

I'm not really clear on how this relates to the dim keyword argument for reductions. Does it? @toivoh, do you have some unifying view of these kinds of dimensions arguments?

[lindahua]

Reduction functions actually do not use keyword arguments. They simply use the second positional argument to specify the dimensions to reduce over.

[nalimilan]

cov is a kind of reduction along the dimension of variables. But indeed dim is ambiguous; vardim looks like a good choice to me, both explicit enough and similar enough to dims used in reductions.

[lindahua]

I will add documentation & news later.