From Statistical Physics to Data-Driven modelling in Biology - session 7

Okay. Um, so we can we can start maybe by recap of what we did this morning before we move on. So we we considered the problem of inferring a graphical uh model which means inferring the graph of interactions between variables and we considered Gausian variables to start with and um well the main results we obtain is that we can do maximum likelihood So if we do maximum likelihood what we have to do is simply solve the equation I mean we simply have to invert let's say the the matrix C the empirical correlation matrix C we invert it we get our estimator for T the true uh matrix that appears in the Gaussian probability and uh and we are we can add a regularization. So if we do let's say if we do maximum likelihood plus L1 regularization then what we have to solve is the equation T minus one uh minus C c minus gamma time the sin of t equal to zero. This cannot be solved explicitly but you can solve it numerically by this algorithm I mentioned called graphical lasso. So it's a very simple iterative algorithm that converges to the solution and um okay uh so we can do maximum likelihood we can do maximum likelihood plus regularization and we can do so likelihood and if we do the likelihood then well I will not rewrite uh maybe all the formula but we can estimate each row of the matrix T So each row which is each row of the matrix is estimated independently by looking at the observed the likelihood of variable I condition to the other variables. And um also in this case we can add uh if we want we can add an L1 regularization. We have seen that for M equal I mean when when the number of data goes to infinity all these estimators converge to the to the correct result. I mean if you for gamma equal to zero if you put regularization then you have something different and and If m is finite then what we have seen is that essentially what you need is what you need is m much bigger than log l and uh a proper choice of regularization gamma and then you can recover the nonzero elements of of t >> then yes >> why can't L1 regular. >> No, no, we can we can >> we can add we can add an L1 regularization. >> No, no, sorry. So, if if you choose if if you have enough data and you choose gamma uh of the order of uh of this then you will recover the non-zero elements of T hat uh with good precision. Okay. Now the next step is what happens if the variables are not gausian. So of course then it depends on what kind of variables you have. But we will we will focus on the simplest uh next step which is uh categorical variables. So discrete variables that take Q states. And this is because also I will uh I will show you an application to proteins and proteins can be proteins are sequences of amino acids. There are 20 amino acids and so you can uh model proteins as categorical sequences with Q equal 20 states or 21 if you also include gaps. So we will discuss this. So now we consider categorical variables. Okay. So categor categorical variables will be discrete variables sigma i uh that take possible states and in this case. So instead of before the var the data before were called y now I call them sigma to distinguish the categorical variables and so we will have each data is a sequence sigma 1 sigma 2 sigma l of length l um okay so if you have um binary data you can choose q equal two but Then typically instead of having sigma i equal to one or two you you you choose 01 or minus one1 it doesn't really matter I mean the name that you give to the states are not very important what is important is that you have two states okay so before we do any graphical model let's see what happens if we have no interaction Okay. So no interaction will be the case where the variables are independent. So if the variables are independent, it means that the probability of the uh L variables is just a product of individual uh probabilities. And uh we will it's always convenient to write these probabilities as exponentials. Uh so we will write this as a product of exponentials of some hi that I will call fields because of the analogy with statistical physics where these terms are usually like in easing models this would be magnetic fields. Um so if I call the probability of sigma I I call it exponential of hi of sigma then I have to normalize so I will divide by the sum over let's say a that goes from one to q of the exponential of h i of a okay so on each on each side I have l sides this is my sequence let's I mean think to a protein sequence so you have each variable represents the state of one one site so on each site I will have in principle I have Q parameters that are these fields okay so sigma i is 1 2 3 so I have one I have h i of of one h i of 2 h i of 3 up to h i of q. Now how do I if I am given data? Now if I have data let's call sigma a collection of sigma let's call capital sigma collection of small sigas. So I have m data points. Each data point is a is a collection of L variables. Okay. So this is a vector. This is in let's say one Q to the L. So this is one data vector and then I have M data vectors that are generated from some distribution. How do I get my H? Well, I I can do once again I can do maximum likelihood. Okay, so let's call this normalization Z I. So this is the denominator. Yes. >> It's exponential of h i of a where a is a is a symbol that I use to sum. So for example, if I have two, let's make an example. If I have two states, okay, in the numerator, I will have for example the probability of one The let's me write it well. The probability that sigma i is equal to one is going to be exponential of h i of 1 divided by z i. And the probability that sigma i is equal to two is exponential of h i of 2 divided by z i. with Z I >> Yeah. Yeah. Exactly. So then now Z I Z I is the sum is what I need to to make the sum equal to one. So Z I is going to be exponential of HI of 1 plus exponential of hi of 2. So I call a an index. So this will be can can you read if I write here? >> Sorry. uh so zed I write it as the sum for a that goes from 1 to 2 of exponential of h i of a which means a is equal to 1 2. So I get a= 1 give me this and a= 2 give gives me this and if I have two states I will have h i of 1 plus h i of 2 plus h i of 3. Okay. So this a means just the state of the variable that I have to sum. Okay. Okay. Now um I want to obtain the h. So what I can do is to maximize the likelihood again. So the likelihood here is going to be 1 / m times the sum of the logarithm of p of sigma m. So I have my data and I sum all the likelihoods of of all the data. So what do I get? Well, the the probability of a data point is a product of probabilities one per each site. So the logarithm is a sum of logarithms. So what I get is the sum over i from 1 to l of the sum over m of the logarithm of p i of sigma i. Okay. So for each site I get the average over the data of the logarithm of pi of sigma i and this what is this? This is one uh this is so I get sum over i. So when I take the logarithm of this is pi when I take the logarithm I get 1 / m times the sum over m of h i of sigma i m. Um when I take the log of this I have h i minus log z i right. Okay. So this is my likelihood. So it is useful to rewrite this in a way that looks a bit more complicated but it is useful which is to write this exponential of hi of sigma i. So I I will write this as a product over I of exponential of sum over A from 1 to Q of HI of A times delta of sigma I A. What does it mean? Divided by Z I. So it what I do here is I sum from one to Q. So I do hi of 1 times delta and this is this is a chronicer delta. So if sigma i is equal to 1 then this will be one when a is equal to 1 and zero otherwise. So I sum h i for example if sigma is equal to one I do h i of 1 * 1 plus h i of 2 * 0. So it's the same as this. Why I want to do this? because then my likelihood I copy what is here. So this h i of sigma i is written in this way. So it will be the sum over i from 1 to l. Then I get 1 / m times the sum over m. And here instead of h i of sigma i put this. So I get sum over a from 1 to q of h i of a time delta of sigma i m a okay minus log z. Why I want to do this? Because now what I get here is sum over i sum over a of hi of a and then I exchange exchange the sums and I get 1 / m sum / m of delta of sigma I m a minus log z what is this what I do here is I I sum over my data and every time for example if a is equal to one I take all my data and every time that sigma i is equal sigma i is equal to one I count one. So this is the total number of ones that I get in that position i divided by the total number of data. So I will call this f i of a. This is the frequency of on site i. This is the frequency of symbol a. A goes from one to to Q. So for each A, basically I have my data. Let's see. This is my this is my data matrix sigma. Each row of the data is one vector of of data points. So if these are proteins, each row is a protein sequence. I take site I and I look in each sequence what is the symbol. So for example here I have one here I have three here I have seven one two six whatever and then I have to count how many times I have one divided by the number of rows and this is just the frequency of having one. Okay. So, so what we get is that now the the likelihood is just the sum over i of the sum over a of h i of a * fi of a minus the logarithm of z I Okay, thank you. So now I have to find the H. How do I find the H? I have to maximize the likelihood. So I have to take the derivative of the likelihood with respect to H I A and set it to zero. What do I get? Well, so first of all, this is a sum over i. So because the sides are independent, each site the h i appear here and here in each term separately. So when I take the derivative, each term any when I derive respect to site one, all the other sites disappear. So the sum over i disappears and then I get here I get fi of a minus the derivative of the logarithm of zed with respect to h i a. What is this? Z i is equal to this this sum. So when I take the derivative I get min -1 / z i times the derivative of z i with respect to h i a which is zed is is a sum of exponentials that is exponential of hy of one plus exponential of h i of 2 and so on. If I take the derivative with respect to hi of two, all the terms except the one with hi of two disappear and what I get is simply exponential of hi of a. Right? What is this quantity? This quantity is the probability according to my model. So my model is this and according to my model this is the probability of observing a in position I. So I will call it pi of a. Okay. And so what you see is that the maximum likelihood equations tell you that the probability according to your model of having amino acid 7 in position I should be equal to the frequency of amino acid 7 in position I that you see in the data. So this is an equation for the hi A. So you have here you have L* Q parameters that are the Q fields on each position and these are L* Q equations that will give you uh these parameters. Um you see that these are moment matching equations because basically the frequency of of I a can be written as the expectation over the data as here. This is the expectation over the data of delta of sigma i a. So this is the this is a function that counts one every time sigma is equal to a. So this gives you the frequency of and this one is the same but in the model. So these are moment matching equations that tell you that the expectation over the model of something should be equal to the expectation over the data of the same thing. And we have seen this already last week when we did maximum entropy. And indeed this model here is an exponential distribution that corresponds to maximizing the entropy under the constraint that the frequencies of the amino acids or of the symbols on a given site should be fixed. Okay, that should be fixed to those in the data. Okay, so this is just to say once again that doing maximum likelihood using an exponential model is the same as doing maximum entropy with the associated observables that in this case are the the the frequencies of of the amino acids. Okay. So, I can erase this just to be a little bit more concrete. >> Yes, sure. >> Yes. Why should it matter which position? >> H because because I'm assuming that the frequencies can be different. >> Well, no, no, these are not probabilities. This is precisely the normalization. So, >> is the probability? >> Well, no. Hi is is log of the probability plus a constant because because the constant is precisely what I'm using to normalize. >> So, the sum is not one. But this is an important point because in fact there is a redundancy here. I I don't have Q independent parameters on each side. I have Q minus one. So there is so one of the fields can be eliminated and for example you could decide that you nor you add the constant to all the fields you shift all the fields in such a way that the Z is equal to one. So I I I was I I will come to this. So you can decide uh how you normalize I mean how you shift the fields um because precisely I mean if you add the constant to all the fields you will add the same constant in the denominator and the constant will will cancel. So uh let let me just see in which order I wanted to say the things. uh yeah I mean we we can discuss this uh now so let's let's consider for just to maybe make it a bit less abstract let's consider the case of binary variables and now if I use binary variables for convenience I will use 01 which in the previous case would be one two but I mean it's just it will only change the sums here that will go from 0 to q minus Um so what happens if you have binary variables that you have two then you have two fields but you can write the probability of sigma i as the exponential of h i of sigma i divided by and then you can uh for example you can remove from the numerator and the denominator hi of zero. So you can write this as so I multiply above and below by exponential of minus hi of zero. So I get one plus right and then I can say okay this term is equal to zero if sigma is equal to zero and otherwise it is equal to hi of 1 minus hi of 0. So if I call simply hi equal to hy of 1 minus h i of zero I can write this as the exponential of h i sigma i divided by 1 + exponential of hi right so this is h i so in the denominator it's okay and in the numerator if sigma is equal to zero I get zero in the numerator If sigma is equal to one, I get this. So you see that out of two fields, you actually only get one field and you can decide how to normalize things. So this is convenient and it's it's one of the possible conventions that people choose. But so this is to show that these models in in categorical models you have this overparameterization that you can use depending on your uh taste or or what you want to do in order to uh write things in different equivalent ways. Okay. So this is what you get when you don't have interactions. Now if we want to make a graphical model we need interactions. So um so how do we write a model with interactions? We do the same we did in gausian variables. So we can add we can make a model where we add a pair interaction. So we write the probability of sigma as exponential of sum over i of h i of uh sigma i plus sum over i smaller than j of some jig j sigma i sigma j. Okay. And then I have to normalize this by a partition function. Why this choice? Well, you can. So this is what is called a Bzma machine. So why this choice? Um you can see this as usual. I mean you you can think that you do uh maximum likelihood. So you can say okay this is my model. I I choose my model for the likelihood of the data sigma given the h and j and so why this model? Well I mean you can motivate it in many ways but for example you can think okay I'm doing I have in mind some something like statistical physics. So I I say this is the exponential of some energy and my energy is the sum of one side terms plus two side terms and then I can add three sides terms and I can keep going or so and then you do maximum likelihood and you get your parameters by maximum likelihood or equivalently you can say I do maximum entropy and I choose some observables and I maximize the entropy. So which one which observables did you did we choose here? So I do the same as before. I write this as exponential of sum over a i a of h i of a time delta of sigma i a. This is the first term. And in the second term I get I do the same trick and I get jig j of a b delta sigma i a delta sigma i. Okay. And then I normalize. So now this is this this is a distribution that is the the exponential of a linear combination where the h i a and the j i j a b are the lranch multipliers and these are the observables you want to fix. So this corresponds to do maximum entropy with the requirement that the expectation over data of delta sigma i a should be equal to the expectation over the model of delta sigma i a. So this is the fields term and the couplings correspond to imposing that expectation over the data of delta sigma i a delta sigma i should be equal to the same over the model. Okay. So what are these things? The first is the frequency or the probability of finding a seven in position I. What is the second? The second is a product of delta. So it will be zero unless sigma i is equal to a. Sorry, there is a j here. This is j. Sorry, because it's sigma j here. So this is zero unless they are both one. If they are both one then it means that sigma i is equal to a and sigma j is equal to b. So for example if I do fig j of three six if I have it means I have to take two sides and I have to count how many times I get three here and six here. Okay. So if there is another uh three here and six. So this will be so this is what is called so this will be the frequencies and these will be something like co occurrences or something like this. So you count how many times you find two symbols together in two positions. Okay. So you can see this choice of model either as maximum entropy model condition to these two u uh constraints or as just a model and then you do maximum likelihood and you get your parameters. Okay. So based on what we said last week uh when you maximize either the entropy if you think in maximum entropy or the likelihood if you think maximum likelihood it's the same if you have no prior I mean it's the same then you you will get the equations for the h and the j that are given by by these conditions okay so how many parameters we have we have l* Q fields and we have for each pair of sides we have a matrix that tells us what is the probability of one symbol here and the other symbol here. So it's a Q square matrix. We here this is going to be this is symmetric if I exchange I with J. But then when I write it in this way, it's not symmetric if I exchange A with B. So for each symmetric pair I J, so I have L * L - 1 / 2, I have Q ^ 2 terms. Okay. Okay. However, these are many parameters but also in this case I have some redundancy because so let's call the frequencies fi of a and let's call this pi of a and if I call this fig j of ab and I call this pig J of AB. What I have is that the sum of the frequencies sum over A of the frequencies should be equal to one. So these are not all independent. There are Q minus one independent constraints and the the the last one will follow from the normalization. And here I have that if I sum this over a, if I sum over a fig j of ab what do I get? I'm asking what is the frequency of getting a free here? If I get anything else here, well then a will disappear here. when I sum over a the delta is always one and what I get is just fj of b. Okay. So also these quantities are not independent. There are two equations like this one for each b and there is another one for for the other side. So there is a sum if I write sum over b of fi j of a b I get fi of a. So I have on each uh site I have one less constraint which gives me the possibility of choosing one field as I want as I did before. And on each pair I have the possibility of choosing Q. So basically one row and one column of the matrix JI J A B. This is a Q by Q matrix for each pair. I can choose it as I want. So I have this freedom of that is called sometimes uh for technical reason is called gauging variance. It's not really a ging variance but it's called like that. Uh so you you you can do basically you can choose a parameterization of your model with this freedom. Okay. >> Yes. >> Explain the freedom again. >> Yes. So the point is that I have Q let's look at the fields. Let's forget the the couplings. I have Q fields on each side and here it looks like I have Q equations one because okay the fields are I look at site at this site and I have the Q fields that are the fields for symbol 1 2 3 Q and then in principle here I have Q equations because I have I have to impose that the probability of the model of having a six here is the same as the probability in the data of having a six here and I have to do it for 1 2 3 up to Q. So it looks like I have Q equations for Q parameters. In reality I have Q minus one equation because FI of A is the frequency of symbol A on site I but the sum of the frequencies must be one because there must be a symbol on site A. So the sum of the FI is equal to one. which means that if I know Q minus one of them, the last one is just one minus the sum of the others. So in reality, these these equations are linearly dependent. One of the equations is a linear combination of the others. So this means I have Q parameters, Q fields and Q minus one equations. So I have the freedom to choose one of the fields. Okay. So I could set to zero one of the fields or I could uh decide as I did before that I call field the difference between two fields. I can play with this and I can choose one of the fields for the couplings is the same. Here I have in principle here I have Q for each choice of I and J. Here in principle I have Q square equations that correspond to all possible choices of the first symbol on site I and the second symbol on site J. But in reality they are not linearly independent because if I sum over a FIG the sum of the delta will be one. There is always one symbol. So I get rid of this and I get the expect and I get the expectation of of this which is the frequency. So the the two point the two point correlations are not Q ^ square in reality they are Q ^ square minus Q minus Q I mean plus one because okay then so basically you have a matrix Q by Q and you can choose one row and one column as you want. For example, you can fix to zero all the elements here. Okay. So you have an over parameterization. This is not I mean it's not that important but in some cases it can lead to we will see we will see in a second that in some cases it's important to have this in mind. Okay. Well, so now how do we get this H and J? Okay. Um, now things are much more complicated because um the problem is that I told you that if you want to do maximum entropy, you can always do Boltzman learning. So you could say okay I I do Boltzman learning. So what do I do? For example, I start from zero. This is my probability of sigma. I start from h= 0, j equal to 0, then it's it's all fine. So I have uniform if I put h equal to 0 and j equal to 0 this is a con. This is exponential of zero which is one. So I gave the same probability to all possible states. I get the uniform distribution. And then I can say okay I have to to solve these equations. So I can write that the derivative of h over time let's say. So I do gradient uh okay maybe sorry I wrote it in a different way the other days. I can say okay my field at time t + one is going to be my field at time t plus a learning rate times the difference between f i a minus p i a. Okay. And I do the same for J. I say okay J I J A B at time t + 1 is equal to jig j at time t plus sorry these are at time t. Uh so you remember maybe from the from last week we said okay the the free energy is uh convex we want to minimize the free energy. So we can take the derivatives and we can do gradient descent. If we do gradient descent, what we get is that at each step of the gradient descent, we update the fields and the couplings. So we have a learning rate and the update is precisely the diff we want to find the fixed point where f is equal to p and fi is equal to pi j. So the update will give us the difference between the two. You can you can take the notes from last week and you can work it out. So this means I start with for example zero j and zero h and then I compute pi of a and pi j of a b in the model. These ones I know them from the data. I I just have to count in my data how many times I get a certain uh symbol in a certain position. So these are computed at the beginning and you don't touch them. At each time step you will update H and J and you then you have to recomp compute PI T and PI JT and keep going. The problem is that computing these two quantities is not easy. Why it's not easy? because in order to compute these quantities you have to know the probability and in particular you have to know the the denominator of the probability and then you have to sum over uh all the other variables. So in order to compute these quantities you should be able to do to normalize these quantity and the partition function is the sum over all the possible sigma and the number of sigma that you have is is uh q to the l right so it's exponential in l it's a huge number you cannot do it so how do you do this in practics and I will show you a code that does it you do Monte Carlo okay so at time t you have your model you do Monte Carlo you generate so these are your data at time t you do Monte Carlo with this uh distribution you generate a similar alignment of data generated from the model so these are your data you generate from your model a similar uh number of uh states from Monte Carlo and then you do the same calculation. You count how many sim how many six you have here and you match with and you try to match with and well you can do that I will show you but it's it's quite heavy because the um the Monte Carlo can be slow you can you can have convergence problems you have to make sure that your data that that the data you generate from the model are well equilibrated and so on and so forth so I don't know if you have ever done Monte Carlo but um if you tried well you know you should see I mean in many cases when you do Monte Carlo you have slow uh convergence okay so okay this is this is doable it's called Boltzman learning but the price to pay is that you have to do Monte Carlo because you cannot compute uh zed okay is it clear up to this point? Yes. >> Yes. Uh so you can see it in two equivalent ways. One is to take the log likelihood. So you say okay this is my model. I take the I do the lo the log likelihood. So if you do the log likelihood you get 1 / m s of log p. So you take the log what do you get? You get uh when you take the log of the numerator, you get sum over i of h i of a times fi of a because this one as before you can write this as sum over i of h i of a delta sigma i a. Okay. So when you take the log and you do the you take the expectation of of this over the data you get this. The other term gives you sum over i j j i j ab f j ab a b for the same reason and then you have minus log z. >> Yes. Thanks. Here it's Thank you. So now when you take a derivative with respect to h i, the derivative of l of h i of a is going to be fi of a here minus the derivative of log z. When I take the derivative of log zed, I get 1 / zed times the derivative of the numerator with respect to h, which gives me the probability. This will be pi of a. Okay. And the same if you do the the derivative with respect to j. This is precisely the derivative of the log likelihood with respect to h. And so this is gradient ascent. You go up the gradient. If you do maximum entropy, you minimize the I mean you maximize the entropy or you minimize the free energy and you get the same result. Okay. Okay. So >> yes. without form assuming that there are no correlations. So why do the same? >> No, it depends on the model. So I give you some data. Okay. and you let's say you don't know where the data come from. So you can say okay let the first attempt is I don't put any correlations. So I I set J equal to zero. In that case I only have the fields and I and the partition function in that case it's easy to compute because it's factorized. So and then you can compute the H and you get your H. This gives you a model but it doesn't give you an interaction. So then the point is if you think that there are interactions you can add them and see what happens but the price to pay is that now you have to do Monte Carlo. What do you do with interaction? So the point is um as I was saying last week I mean when you do inference you have to you need to have a model in mind of what you want to infer. So if you do this you you can see it in two ways. One way is to say okay I want to do maximum likelihood. I have to choose a likelihood. This is the minimal likelihood that has one point um interactions and two point interactions. So this is a graphical model like in the gausian case where I have terms that connect two variables. So if I write the conditional probability of site I, it will be the exponential of the field plus a sum of individual couplings of I with J. I could add three point interactions for point interactions then it's up to you where you want to stop. The more you add things and the more you have parameters and the more it's difficult to fit your model. So then it's a choice. Otherwise equivalently you can think okay I don't know what to do. I don't know my model but I want to maximize the entropy of my model condition to the data. What what do I want to see in the data? If I'm looking for a graphical model, I can say okay what I want to what I care about is individual uh sites and pair interactions. So I will I will do maximum entropy condition to the fact that the model should reproduce the single site statistics and the two side statistics and then you get the same. So you can see it in these two ways that are equivalent and of course you are throwing away something because there could be free body interactions and then the problem is that the free body interactions are going to be L cq times Q cq. So if you have a protein with 20 amino acids and 100 sides you have a big number of parameters if you want to include all of them. So in practics in in all the applications I've seen if you try to to do more than this you will typically have to make strong assumptions on the free body interactions that maybe they are very sparse or uh or you put a strong regularization or things like that. So this is a minimal let's say the model without the couplings is the minimal model that you can do simpler than that is doesn't make much sense and this is the minimum you can do to introduce a correlation between variables let's say and so again the aim here is the same as in the case of gausian variables we want to infer the J to infer a network of interactions between sides and these J are going to be the direct interactions that are um and we will see what kind of information we can extract from them. >> So we can also >> yes so we we can add regularization and we can do exactly the same kind of things. So we we I will come to that. Okay. Um so up to this point is it okay? So this is let's say this is a simple framework where you do a simple uh let's say set of assumptions for example you say okay I do maximum entropy condition to the single site and two sight statistics I get this Bzman distribution this is also called Bza machine I want to get the parameters I run my Bzman learning I get the parameters and then I can see what kind of information I can extract but I will have some network of jig that gives that gives me the interactions between these variables. The problem as I was saying is that this is a bit um slow. I mean I will show you a code that can do this for proteins in maybe one hour on a GPU. So you now you can do it on your uh computer but uh if you go to big systems and so on you can have problems. So there are alternative u schemes that are simpler. So I would like to before I go into the problem of regularization. So for the moment let's say we don't put regularization. So before I go into the problem of regularization I would like to discuss alternative schemes to obtain this J and H more efficiently. Okay. So to do that however let me I hope this will not confuse you but let me restrict to a simpler case. So let's take again binary variables. Okay. So we go again to binary variables. You can do it also for Qstate variables but it's a little bit uh just the notation is a bit heavier. So I put binary variables with sigma equal to 01. How does this model look like? So my P of sigma. So as I was saying before, I can write this as I can say okay I have my HI of sigma which means I have hi of zero HI of one but I have this freedom. I can choose one of the fields. So I will choose this to be zero or or if you want I I will substract subtract this from both and I will write it in in this way zero for state zero and hi. So this means here I get exponential of sum / i of h i sigma i. Now what we are going to do? So we are going to see first I will specialize these two binary variables and then I want to show you alternative ways to get the h and the j without doing the boltsman learning that is computationally heavy. application to icing. >> Yes. So I mean well yes this is the the icing case. Yes. When they are binary and you can think that for example the neuron data that I showed you last week we we we binarized it. So we could use you could use this to describe neurons when you have two states. So basically are we expecting that to learn the magnetization >> the magnetization function? >> Um well this is more than the magnetization here we want to learn all the local fields and to body interactions. So okay let me write and then I I so okay let's specialize these two binary state two binary variables. So I have this term will become if I if I assume that I have zero field for zero and hi for one then I can write it in this way because if sigma is equal to zero I get zero and if sigma is equal to one I get hi and similarly for the matrix J. The matrix Jig J is a matrix J of 0 0 J 01 J 1 0 J11 1. But I told you that I can choose one line and one column. So I will put everything to zero except this. I will write it as 0 0 Jig. Okay. And then I get sum over i j of jig j sigma i sigma j. So this is my model and this is an ising model with local fields hi but these are I have a diff I'm assuming that I have a different local field on each site and then I have couplings jig j that are different for each pair and I would like to get the h and the j. Okay. So why for why why why do I want to do that? Just to make an a concrete example, you remember the data I showed you of the neural recording in the rat uh during um sleep and then uh task and then sleep. No. And the goal was to infer net a group of co-activating neurons that were um co-activated during the task and then were replayed during the sleep. Okay. So you can think that your data is well okay the the the neuron data were transposed with respect to the but you can okay this is time okay and these are the neurons so you you had 37 neurons So each line here is one is the evolution of a neuron during time and the state at each time can be zero or one. You can say zero is the neuron is silent and one is when the neuron is um active. So what you want to do you can say that at each time this is your sigma vector okay and so at each time you have one realization of sigma and you can say okay let's assume that they are all independent realizations so these are my data so I say this is sigma 1 this is sigma 2 sigma 3 sigma M. So these are each line is a vector of data and each line contains 37 spins. So I I can try to fit this model to this data. So I I I maximize the likelihood on this data. I do BSA learning. Okay. And then I will get for each neuron I get a field hi and for each pair of neurons I will get a coupling jig and the hi will tell you how active this neuron is. So the bigger hi the more the neuron is active. So we have seen in the data that there were neurons that were very active and neurons that were very silent. So this is an individual property of each neuron while the jig tell you what is the likelihood of the two the two neurons being active together. The the bigger is JJ and the more likely it is that sigma and sigma j are both one. Okay. And so you can apply this method to to to those data and you will find that the neurons that were um identified by the PCA are also those that have strong JIG they are strongly coactivated. So there is a paper by uh by the group of Remy Monos and Simonoko and I think Ferrari where they did where they analyze those data using this method. Okay. So the goal is it's so it is an anime model but the the usual setting of the model in physics is that the field is the same for everyone and the J is nearest neighbor formanatic coupling and then you want to compute from this you want to compute correlation. Here the setting is that from data you want to infer this and because the data are uh I mean because the for example the the spins represent different neurons they can have different fields and different couplings and so it's it's all going to be disordered. Okay. I I hope it's clear but feel free to keep asking questions. uh uh as much as you want. Okay. So now the goal is how do I get I'm given the data. How do I get H and J? So one way is Bzman learning. Boltsman learning is exact. It will converge because the function is convex. So it if you wait long enough it will converge. But it is slow because you have to do Monte Carlo to sample from this model and you have to repeat it many times. So we want to see if we can do something faster and uh and cheaper. So we will do approximations. Okay. So the first approximation is what is called minfield. So what what do you do in mfield? Okay. So what you do in minfield there are many ways of doing minfield. Uh but one simple way is the following. So you start from an exact relation which is the following. If you compute the magnetization of spin I. So the average of sigma i now this this average is the average over this distribution. So this is the sum over sigma of p of sigma times sigma i. Okay. But you can also write this as the sum over sigma of p of sigma i condition to the other spins times the probability of the other spins times sigma i. Right? What is the probability of sigma i condition to the others? We will need it. So let's write it now. Now I do exactly the same thing as I did for the gausian case. What is the probability of sigma i condition to the others? I have to I extract from here all the terms that involve sigma i and the others are constants. Right? So this will be proportional to what I have exponential of h i sigma i from here and then the other fields are not important because the other sigma the other sigma are constants so are fixed I'm conditioning on them plus here I have what I have to sum over i and j on all possible pairs but I only want the pairs that involve i. So this will become a sum over j different from i of j i j sigma i sigma j and all the pairs that do not contain i. So all the j k where there is no i are constants. So this is my probability that I have to normalize. How do I normalize it? I have to divide by the sum of this over sigma i. What is the sum over this over sigma i? Because sigma is binary. This is I have to sum the numerator when sigma is equal to zero and when sigma is equal to one. When sigma is equal to zero, I get zero. So I get the exponential of zero which is one. Let me let's simplify this. I can factoriize here a sigma i. So I will get this times sigma i. So when sigma is equal to zero I get one and when sigma i is non zero is is one I get I get this right. So this is the conditional probability and I can write it as the exponential of a field hi time sigma i divided by 1 + hi where hi is this this term in parenthesis. Okay. So now let's go back here. The average of sigma i over the boltsman distribution is this. I can write it as the sum over sigma. Sigma means all the sigma of p of sigma i condition to the others times the probability of the others times sigma i. Now this can be written then I can write as the sum over all the other spins but i of this And then I have the sum over sigma i of p of sigma i condition to the others times sigma i. What is this? It is the average of sigma i conditioned to the others. So I can compute it from here. So if I multiply this by sigma i and I sum the denominator is this and the numerator is the sum of sigma i times this which is just exponential of sigma. So I get sum exponential of hi. So I get this times exponential of hi divided by 1 + exponential of hi. Okay, this is the this average here is given by this where hi is this thing is the h the small h i plus the sum of all the interactions with all the other spins. So this is the average of h of this thing. Okay. Where where now the average is over this but actually hi does not depend on sigma i. So taking the average over the probability of the sigma minus i or the full probability is the same because this does not depend on sigma i. Okay. So I get an identity. This is an exact identity which tells tells me that the average of sigma i is equal to the average of the exponential of h divided by 1 plus the exponential of h where h i is the the field plus the sum of the interactions with the other spins. Okay. Is this clear up to this point? >> Yes. >> Which one? So, okay. I I I raised something, but up to this point, it's it's okay. So, I I I say, okay, I have to average sigma. So, I write the probability of all the spins as the product of the probability of the ones that are not I times the probability of sigmi condition to that time sigma i. Now I do this average. This average is this is the probability of sigma i condition to the other spins and it has this form. When I multiply, if I multiply by sigma i and I sum, in the denominator I get 1 plus a this and in the numerator I get if if sigma is equal to zero, I get zero. And if sigma is equal to one, I get e to the h. So I get this. Okay, sorry, sometimes I skip some steps, but so I get this. And now I say okay now this has to be averaged over the probability of all the other spins except I. So this is the marginal probability where I sum over sigma i. But because this thing does not depend on sigma, it doesn't change anything. I can also say this is the I can replace this by the sum over all sigma because this is not this is independent of sigma. So if I write it in this way when I sum over sigma I get the marginal probability of the others. Okay. So I get this identity. Anyways, this is an exact identity called these are called column identities and well these are exact identities. Okay. So I can I erase? Okay. So the the mean field approximation can be obtained in many ways and one way is to to take this and just approximate this this average of this complicated thing. We can approximate by replacing this by 1 + the exponential of the average. Okay, so it's an approximation. How do you justify this approximation? You know it. I guess the idea is that hi is given by a sum of many terms. I mean it's given by a sum of terms and the number of terms is given is given by the connectivity of I. So if if if this is spin I, every time spin I is connected to another spin J, there will be one term in the sum. So here you have a sum of terms that are and the number of terms that you have is the connectivity of spin i in the graph that we are trying to infer. If the connectivity is large enough then by central limit theorem this thing will have small fluctuations and so you can replace the average of this function by the function of the average because the fluctuations are small. So the this approximation is good enough when you have many terms in this sum. So for example, if you are in a system in very high dimension where you have many neighbors or if you are on a graph that has large connectivity for example a fully connected graph or a graph that is close to be connect well connected okay otherwise it's just an approximation. Okay so then what do I get here? I get I get let's say well okay the average of hi that appears here is given by h i plus the sum over j of j i j times the average of sigma j I take the average and this is the random term okay so now I can at from this a simple equation for J by the following observation that the derivative let's take the derivative of sigma i with respect to h k. Okay. So if I take the derivative of sigma i with respect to h k. So what I'm what I'm trying to do is I look I change the field on one spin and I look to how the magnetization of another spin is changing. I can take the derivative here with respect to hk and I can take the derivative here. So let's call this function f of the average of h. So f is the sigmoid but I don't need it actually. So when I take the derivative with respect to h k here um I take the derivative with respect to h kk what I get is the derivative of f with respect to its argument. So I get frime of h i times the derivative of h i with respect to h kk. So this is given by frime of hi. When I take the derivative of h i with respect to h k small what do I get? I get the derivative of small hi with respect to small hk which is a delta function. It's delta i k plus sum / j of jig j times so let me call this r i k. So this is the response of spin i to field k. Here I get the derivative of sigma j with respect to h k. So I get R J K right. And so in the end what we get is that this matrix R the matrix R is equal to F prime of let's let's call this DI this term I call it DI. So this is d i times delta i k plus the sum / j of j i j r j k. Okay. So this means that the matrix R is equal to this matrix D. This is a a diagonal matrix plus D times J * R. Okay. And so I can get J from this. So what is J? I multiply uh everything by D. So I will get d minus one * r - 1 is going to be j * r and so then I multiply by r -1 and I get d -1 - r -1 = j Okay. So in the min field approximation the J matrix is equal to uh D minus one minus R minus one here. Okay. Here J is different from I. I dropped it but I assume that J I I is zero. So the diagonal terms are zero. So this means that uh the matrix element jig is equal to so if I is different from j the matrix element jig j this is a diagonal matrix so the of diagonal there is nothing this is the this diagonal so it's simply equal to minus the inverse of R. Okay, what is R? R is the response function. So R is the matrix of the derivatives of sigma I with respect to HK. But so if I assume that this is the probability of my variables, the derivative so can I can I erase this? the derivative of sigma i with respect to h k. You can do the derivative. You take this, you write the average of sigma i. When you do the derivative, another sigma drops and this is sigma i sigma k minus sigma i sigma k. So this is called the fluctuation uh dissipation relation. It's relation between response and correlation. It's the static version of the fluctuation distribution relation. So this is CI. It's the matrix of correlations that you have in your um system. And so what we found in the end from the mean field is that the matrix J that is the matrix we want to infer is just the inverse of the matrix of the correlations. So same as gausian. Okay. In the gausian case we call it t and there was a minus sign. This is just conventions. In the spin literature usually you put the plus here. In the gausian literature usually you put the minus. So J is morally the equivalent of minus T and in the Gausian case we found that T was equal to C minus one and here we find that J is equal to minus C minus one. Okay. >> Yes. >> How did you say that? >> Yes. So for this you have to do a calculation that maybe I will not do but so you write that sigma i is equal to the sum over sigma of sigma i times the exponential of uh this divided by z. And if you take a derivative of this with respect to sigma k to h k what happens is that here there is the sum let's let's call this index k just to avoid confusion when I take the derivative with respect to h kk there is a sigma k here so I get sigma sigma k which is this and then but then there is h k also in the denominator So you get another term which is minus 1 / z² times the derivative of zed with respect to h k times this and with with a little bit of patience you can check that this is precisely sigma sigma k okay I it's better if you do it by yourself because it's very simple but if I write it with all the indices it will be it's this kind of things if you do it once then you're happy and you >> calculation you assume that matrix J is independent of H, right? >> Yes, because my parameters are J and H and I want to learn both. So here I I forgot I I mean the the the H will come from the the in field the the H will come by fixing the the the first moments like if if there was no J. >> Yeah, you said that we should assume that are small No no no not the sorry in which step no this H no which >> small >> no the fluctuations are small when I did the mean field approximation. So I have an exact relation that tells me that the average of sigma i is equal to the average of this. Okay, this is exact but I cannot use it. What I want to do is to approximate this by the by replacing H by its average and I can do this if the fluctuations of H are small not H the fluctuations. H is a random variable because H is this sum and the sigas are random variables. So h is a random variable. If the fluctuations of h are small then I can replace I can approximate a function of h. Sorry. So if the fluctuations are small I can say okay here I have the average of a function of h and I can approximate this by the function of the average. So we don't need that small. >> No no no no no no. What we need to assume is that the fluctuations of this are small which is true if you have many many terms in this sum typically. So if you have a graph that is that has big connectivity then this is a reasonable assumption but it's an approximation. So it's not guaranteed to be exact in any case. The thing is that it's an approximation, but it's very fast because now you don't have to do Monte Carlo BMA learning. They use the GPU blah blah blah. You just do the inverse of a matrix and you're done and you get your J. Okay. >> Sorry. >> Yes, sure. >> Yes. Yes, you're right. So, but what we are going to do then is to approximate C with the data. J, but I want to fit my model to the data. So I have to ch what I want to do is to choose J in such a way that this these things are equal to the data right. So I want this to be equal to the data and I want this to be equal to the data. So I will assume I will say that the C here is has to be chosen in such a way that C is the one of the data and then J is the inverse of this matrix. Okay. So if you want I have I have one equation says that J should be equal to C of J and I have another equation that tells me that C of J should be equal to C of the data. So I replace the second equation here and does it make sense? Okay. Okay. Um I think we can stop here and what I will do to so the last thing I want to do tomorrow before we move to the application it will be more more relaxed with no equations anymore. I want to show you the likelihood calculation. So we can do in this case the same thing we did for gausian variables. So we can take this conditional probability and we can use the conditional probability to learn the J and it will have some advantages over the Boltzman learning. So tomorrow we compare the likelihood with the Bzman learning and then I want to quickly discuss the problem of regularization and then we can finish the theory and uh then the the rest of the time I will show you some uh some application to the data.