Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Discussion] Design of Random Matrices #17039

Open
czgdp1807 opened this issue Jun 16, 2019 · 17 comments

Comments

Projects
None yet
2 participants
@czgdp1807
Copy link
Member

commented Jun 16, 2019

References

Documentations and References
[1.1] https://en.wikipedia.org/wiki/Random_matrix
[2.1] https://arxiv.org/pdf/1712.07903.pdf
[3.1] https://www.wolfram.com/language/11/random-matrices/

Associated PRs
[1.2] #16866

Aim

This issue is raised for the following purposes(in the order as they are mentioned):

  • Discussion of API for Random Matrices.

  • Discussion of the class structure for Random Matrices and its associated features.

We will first dicuss API and after completion we will move to the discussion of class structure, if needed.

Checkpoints

[1] Use Cases for API - #17039 (comment)
[2] Proposed API - #17039 (comment)
[3] Final API - #17039 (comment)
ping @Upabjojr

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jun 16, 2019

The use cases for API are mostly taken from 2.1. For familiarity of the reader, random matrices are those matrices which have at least one random variable in its entries. The intext problems in 2.1 suggest the following:

Note - H is a random matrix with its entries as iid random variables.

  1. What is the jpdf ρ[H] of the N**2 entries {H[1,1] , . . . , H[N,N]} of the matrix H in (1.1)? - The answer is simple as all the entries iid, its the product of the independent densities.

  2. What is the jpdf of the N(N − 1)/2 entries in the upper triangle of the symmetric matrix Hs in (1.2)? - Again the answer is simple, refer, the page 10 of 111 or may be 11 in document viewer. Hs is (H + H_transpose)/2.

  3. What does the jpdf of eigenvalues ρ(x1 , . . . , xN ) of a random matrix ensemble look like? - The joint distribution of the eigen values is very much important in random matrix theory. Infact, spacing distributions are also important.

Apart from the above questions, the Gaussian ensembles are also to be implemented. Their theory is given here https://en.wikipedia.org/wiki/Random_matrix#Gaussian_ensembles.

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jun 16, 2019

I propose the following API,

  • First point is to use MatrixSymbol for Gaussian ensembles as they will easily allow defining densities. We can keep the API for the ensembles very similar to other distributions as they are analogous to distributions of random matrices.
>>> R = GaussianOrthogonal('R')
>>> density(R)
# the density of GOEs.
>>> R.joint_eigen_distribution
# the JointDistributionHandmade object
  • For matrices with RandomSymbol as their entries, we can use the current Matrix API with some extra methods for supporting the joint distribution queries.
>>> M = Matrix([[N1, N2], [N3, N4]]) # Ni is normal iid.
>>> M.joint_distribution(M[0, 0], M[0, 1])
# the required joint distribution of N1 and N2
>>> M.joint_eigen_distribution()
# the distribution of eigen values
@Upabjojr

This comment has been minimized.

Copy link
Contributor

commented Jun 18, 2019

First point is to use MatrixSymbol for Gaussian ensembles as they will easily allow defining densities. We can keep the API for the ensembles very similar to other distributions as they are analogous to distributions of random matrices.

You have already defined an IndexedRandomSymbol. Can we use that same object for matrix operations?

The difference between indexed symbols and matrix symbols is that indexed symbols require an index, and matrix multiplications have to be represented with a Summation operator of a summed-over symbol. Matrix symbols on the other hand are implicitly indexed, i.e. they don't require you to access them with indices.

For matrices with RandomSymbol as their entries, we can use the current Matrix API with some extra methods for supporting the joint distribution queries.

Maybe we should define the methods as external. I mean, instead of M.joint_eigen_distribution(), use joint_eigen_distribution(M). In this way we don't need to subclass the matrix class.

@Upabjojr

This comment has been minimized.

Copy link
Contributor

commented Jun 18, 2019

Also consider using sympy.multipledispatch to dispatch on types, it's vaguely similar to operation operloading:

@dispatch(MatrixBase)
def joint_eigen_distribution(x):
    ...

@dispatch(RandomSymbol)
def joint_eigen_distribution(x):
    ...

You can define the same function dispatched on the type of its input parameters.

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jun 18, 2019

You have already defined an IndexedRandomSymbol. Can we use that same object for matrix operations?

Yes, we can use that but with some changes in it to support various features of MatrixSymbol as there is a bit of difference between Indexed and MatrixSymbol. I will try if we can work it out, if not then I will go with MatrixSymbol.

Also consider using sympy.multipledispatch to dispatch on types, it's vaguely similar to operation operloading:

I will have to explore sympy.multipledispatch. Will get back to you on this soon.

@Upabjojr

This comment has been minimized.

Copy link
Contributor

commented Jun 18, 2019

Yes, we can use that but with some changes in it to support various features of MatrixSymbol as there is a bit of difference between Indexed and MatrixSymbol.

Yeah, maybe we need to subclass MatrixSymbol. I would try to avoid subclassing other stuff (like MatAdd, MatMul), they should remain they ones defined in the matrix module.

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jun 19, 2019

I would try to avoid subclassing other stuff (like MatAdd, MatMul)

I don't think MatrixSymbol inherits MatAdd, MatMul we can safely use it. Am I right?

@Upabjojr

This comment has been minimized.

Copy link
Contributor

commented Jun 19, 2019

I don't think MatrixSymbol inherits MatAdd, MatMul we can safely use it. Am I right?

Right, it does not.

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jun 24, 2019

Also consider using sympy.multipledispatch to dispatch on types, it's vaguely similar to operation overloading:

I went through this and I think it's quite nice idea. Infact, As far as I can figure out, it can avoid subclassing MatrixBase. Am I right?

Yeah, maybe we need to subclass MatrixSymbol.

Experimenting on this.

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jun 24, 2019

I experimented with subclassing MatrixSymbol and found out that, MatrixExpr whose method call, MatMul, MatAdd affecting the stochastic_process_types.py. So, I thought of changing the MRO at runtime using the strategy given here. I don't know how acceptable it is, at least changing the MRO isn't acceptable to me. The code is given below,

class MetaClass(type(RandomSymbol)):
    """
    Used for resolving the undesired MRO
    for the class RandomIndexedSymbol.
    """
    def mro(cls):
        return [cls, RandomSymbol, MatrixSymbol, Expr, Basic,
                EvalfMixin, object]

class RandomIndexedSymbol(RandomSymbol, MatrixSymbol, metaclass=MetaClass):

    def __new__(cls, idx_obj, n=None, m=None, pspace=None):
        if isinstance(idx_obj, Indexed):
            return Basic.__new__(cls, idx_obj, pspace)
        if isinstance(idx_obj, string_types):
            idx_obj = Symbol(idx_obj)
            if n == None or m == None:
                raise ValueError("Shape of the random matrix is missing.")
            n, m = _sympify(n), _sympify(m)
            obj = Basic.__new__(cls, idx_obj, n, m, pspace)
            return obj
        raise TypeError("Invalid parameters passed, "
                "expected an indexed object or name and shape.")

    symbol = property(lambda self: self.args[0])
    name = property(lambda self: str(self.args[0]))
    key = property(lambda self: self.symbol.args[1])

However, problems are still there. I would say we should develop Random Matrices by using MatrixSymbol without involving RandomIndexedSymbol. If I misinterpreted your idea of sub-classing then please let me know. Thanks.

@Upabjojr

This comment has been minimized.

Copy link
Contributor

commented Jun 24, 2019

Omg, please don't try to hack the MRO, it's really dangerous!

RandomIndexedSymbol should not inherit MatrixSymbol. They are separate things, RandomIndexedSymbol is supposed to be a scalar, because A[i, j] is a scalar (even if i and j are symbolic).

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jun 25, 2019

Omg, please don't try to hack the MRO, it's really dangerous!

I understand your point. I was just trying to make it work, but it doesn't, __mro__ never changes probably. Well learnt something new. I will not try to manipulate MRO.

RandomIndexedSymbol should not inherit MatrixSymbol. They are separate things, RandomIndexedSymbol is supposed to be a scalar, because A[i, j] is a scalar (even if i and j are symbolic).

Ah! I misinterpreted your point of sub-classing. So, may be you suggested to make a new class, something like, RandomMatrices which will inherit MatrixSymbol. If that's the case then it is easy, I will
do that and will update the discussion with a new check point, with some concrete ideas.

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jul 3, 2019

ping @Upabjojr
I thought over the idea of using sympy.multipledispatch and inheriting the MatrixSymbol in RandomMatrix. Some concrete ideas on API and class design are given below,
API For Gaussian Ensembles

>>> R = GaussianOrthogonal('R')
>>> H = MatrixSymbol('H')
>>> density(R)(H)
# density of GOE
>>> joint_eigen_distribution(R) # launched using sympy.multipledispatch
# the JointDistributionHandmade object

API for Matrices with RandomSymbol as their elements

>>> M = Matrix([[N1, N2], [N3, N4]]) # Ni is normal iid.
>>> joint_eigen_distribution(M) # launched using sympy.multipledispatch
# the distribution of eigen values
>>> spacing_eigen_distribution(M)   # launched using sympy.multipledispatch
# spacing distribution
>>> joint_pdf(M) # launched using sympy.multipledispatch
# joint pdf of all the entries

I have removed, M.joint_distribution because it doesn't holds any meaning. Anyone can create custom joint distribution using, JointDistributionHandmade by passing the pdf, instead joint_pdf(M) will be better to compute the PDF. So, just for two-three methods creating a new class is not good according to me. Instead, external functions, like, spacing_eigen_distribution, joint_eigen_distribution are better to use.

Class Design for Gaussian Ensembles
I think that, we can implement a base class that can be inherited by all the ensembles. We can add methods like density, etc to it. Now, for compatibility, we may need a PSPace for the ensembles, something like, RandomMatrixEnsemblePSPace which will allow calls like, density(R)(H) and will transfer the queries to specific classes.

@Upabjojr

This comment has been minimized.

Copy link
Contributor

commented Jul 3, 2019

I thought over the idea of using sympy.multipledispatch and inheriting the MatrixSymbol in RandomMatrix.

Multiple dispatching is a way to define a method over:

  1. classes that are not related by inheritance relations.
  2. dispatching on more than one arguments.

If there is an inheritance and just one argument, one should always try to use ordinary methods first.

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jul 4, 2019

If there is an inheritance and just one argument, one should always try to use ordinary methods first.

Ah! right, spacing_eigen_distribution and joint_pdf don't need sympy.multipledispatch.
Other than that, are there any more suggestions for the API?

@Upabjojr

This comment has been minimized.

Copy link
Contributor

commented Jul 5, 2019

Why is there a difference between joint_distribution and joint_pdf?

@czgdp1807

This comment has been minimized.

Copy link
Member Author

commented Jul 6, 2019

Why is there a difference between joint_distribution and joint_pdf?

I thought that returning pdf will be better as in this case it will be more useful, though including both will also be good. One can call anyone of joint_distribution, joint_pdf according to the needs. Well, I think that's not an issue, we can decide on it later on.
Anything other than this point?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.