Reason why f(X).mean and f_post(X).mean do not have the same datatypes?

[patel-zeel]

I was trying out the following code

from stheno import GP, EQ     
import numpy as np

X = np.random.rand(10,1)
y = np.random.rand(10,1)

f = GP(EQ())
f_post = f | (f(X), y)

print(type(f(X).mean))
print(type(f_post(X).mean))
# print(type(f_post(X).mean.mat)) # This returns <class 'numpy.ndarray'>

output

<class 'numpy.ndarray'>
<class 'matrix.matrix.Dense'>

I was wondering that what could be the benefit of returning different data types at cost of consistency?

[wesselb]

Hey @patel-zeel,

You're right that the wrapping in another type is not necessary here.

The reason why these types pop up is to allow structured representations of matrices. For example, consider

from stheno import Measure, GP
import numpy as np

with Measure():
    f = GP(1) + (lambda x: x) * GP(1)

x = np.linspace(0, 100, 100_000)

Then

>>> f(x).var
<low-rank matrix: shape=100000x100000, dtype=float64, rank=2>

>>> import lab as B

>>> B.diag(f(x).var)
array([1.000000e+00, 1.000001e+00, 1.000004e+00, ..., 1.000060e+04,
       1.000080e+04, 1.000100e+04])

which would not be possible if structure was not preserved. Moreover, the type Dense behaves just like a normal NumPy array: it knows how to matrix multiply. The catch is that this does require the use of B (B.matmul, B.squeeze, et cetera), which is suboptimal.

Perhaps the variance is allowed to be a structured matrix type, but the mean should always be a normal tensor? Would that sound reasonable?

[nipunbatra]

HI @wesselb,

Amazing work on Stheno and related tools! The code is perhaps the most natural I’ve seen in any GP framework.

@patel-zeel and I were having a discussion on the above issue. I was primarily approaching Stheno as a user. I wrote the following code to plot the GPs (both prior and posterior).

def plot(f, x):

    mu =f(x).mean.reshape(-1, 1).flatten()
    var = np.diagonal(f(x).var.mat)
    std = np.sqrt(var)

    plt.plot(x, mu, color='purple', lw=4)
    plt.fill_between(x, mu-std, mu, alpha=0.1,color='purple')
    plt.fill_between(x, mu+std, mu, alpha=0.1, color='purple')

    plt.fill_between(x, mu-2*std, mu, alpha=0.1,color='purple')
    plt.fill_between(x, mu+2*std, mu, alpha=0.1, color='purple')

However, as you’d know this will work for the prior and not the posterior since the mu in posterior is accessible as: f_post(x).mean.mat

I was suggesting if the same interface can be provided to access both the prior and the posterior mean and variance.

[wesselb]

Hey @nipunbatra,

Thank you for your kind words. :)

Firstly, I should say that using *.mat is not a safe way to get the underlying tensor. The recommended way of doing this, which will always work, is to use B.dense: B.dense(f(x).var) rather than f(x).var.mat. (For example, if f(x).var is of type matrix.LowRank, then f(x).var.mat will break.) Perhaps this should be emphasised in the README.

Secondly, for predictions, the above code can simplified by using Normal.marginals or Normal.marginal_credible_bounds:

mu, var = f(x).marginals()

or

mu, lower, upper = f(x).marginal_credible_bounds()  # lower = mu - 1.96 * std, upper = mu + 1.96 * std

Thirdly, if you wanted to do it manually and you happen to already be using LAB (the linear algebra inferface behind Stheno), then you should be able to write the above code without any worry:

import lab as B

def plot(f, x):
    mu = B.squeeze(f(x).mean)
    var = B.sqrt(B.diag(f(x).var))

    plt.plot(x, mu, color='purple', lw=4)
    plt.fill_between(x, mu-std, mu, alpha=0.1,color='purple')
    plt.fill_between(x, mu+std, mu, alpha=0.1, color='purple')

    plt.fill_between(x, mu-2*std, mu, alpha=0.1,color='purple')
    plt.fill_between(x, mu+2*std, mu, alpha=0.1, color='purple')

Taking a step back, perhaps it is more natural to let matrix.AbstractMatrix implement __array__ (see this) so that those types are immediately compatible with NumPy functions. In that case, we wouldn't have to worry about converting anything and use f(x).mean like it was a NumPy array. Perhaps that's the best solution here?

[wesselb]

I've added a few methods to matrix.AbstractMatrix which makes the following code run as expected:

import numpy as np
import matplotlib.pyplot as plt

from stheno import GP, EQ

def plot(f, x):
    fdd = f(x)
    mu = fdd.mean.flatten()
    var = fdd.var.diagonal()
    std = np.sqrt(var)

    plt.plot(x, mu, color='purple', lw=4)
    plt.fill_between(x, mu-std, mu, alpha=0.1,color='purple')
    plt.fill_between(x, mu+std, mu, alpha=0.1, color='purple')

    plt.fill_between(x, mu-2*std, mu, alpha=0.1,color='purple')
    plt.fill_between(x, mu+2*std, mu, alpha=0.1, color='purple')

    plt.show()


f = GP(EQ())

x = np.random.randn(10)
y = f(x).sample()

f_post = f | (f(x), y)

plot(f_post, np.linspace(0, 10, 100))

Would this be satisfactory?

[nipunbatra]

This looks good. Thanks for the quick response @wesselb

Thanks for clarifying about .mat. I think I used it as I’m accustomed to looking at available methods in a class looking at the dot completes. I agree it would be a good idea to clarify in the read me about this.

The example looks fairly neat now!

[wesselb]

Great to hear that it looks fairly neat now. :)

I've extended the 20 seconds example to make this point more obvious.