# Optimizing stan code runtime for matrix factorization

**URL:** <https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492>\
**Category:** Modeling\
**Created:** [March 14, 2024, 8:07pm UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492 "2024-03-14T20:07:26Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Garren\_Hermanus](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/garren_hermanus/32/17356_2.png) [@Garren\_Hermanus](https://discourse.mc-stan.org/u/Garren_Hermanus)\
**Post date:** [March 14, 2024, 8:07pm UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492/1 "2024-03-14T20:07:26Z")

</div>

Hi all

I am currently working on implementing some tensor/matrix completion framework but by factorizing the (square) matrix slices from the tensor into lower rank matrices and implementing some correlating structure between the matrix slices. That is let the tensor of partially observed entries be denoted by \boldsymbol{\mathcal{A}} \in \mathbb{R}^{T \times M \times N \times N}. We can factorize the N \times N matrix slices of \boldsymbol{\mathcal{A}} into a collection of lower rank matrices, \boldsymbol{\mathcal{U}} \in \mathbb{R}^{T \times M \times D \times N}, \boldsymbol{\mathcal{V}} \in \mathbb{R}^{T \times M \times D \times N} , arranged in a tensor form, such that:

\boldsymbol{A}\_{t,m,:.:}= \boldsymbol{U}\_{t,m,:,:}^{T} \boldsymbol{V}\_{t,m,:.:} + \varepsilon

We can look at the dimensions t, m as being accompanied with some input data, lets denote this by \boldsymbol{x}\_{a,T}, \boldsymbol{x}\_{a,M}. We want to introduce dependencies along the entries at positions (row and column) ij and ji across all the matrices by using some GP model. That is we denote the vector of entries extracted across all t, m at positions ij, ji as \boldsymbol{a}\_{ij} and \boldsymbol{a}\_{ji}. We combine these into a single vector such that \boldsymbol{a}^{(ij)}=[\boldsymbol{a}\_{ij}, \boldsymbol{a}\_{ji}].

Now for the modelling part. Suppose we have some priors on \boldsymbol{\mathcal{U}} and \boldsymbol{\mathcal{V}}, P(\boldsymbol{\mathcal{U}}), P(\boldsymbol{\mathcal{V}}). And we derived some prior for structure which incorporates the belief that the entries in \boldsymbol{\mathcal{A}} is from a GP and from the matrix factorization, that is prior densities specified as P(\boldsymbol{a}^{(ij)}|\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}}). The probability density of all the reconstructed entries can be given as the product of the individual ij's, that is:

P(\boldsymbol{\mathcal{A}}|\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}}) = \prod\_{i=1}^{N} {\prod\_{j=i+1}^{N} {P(\boldsymbol{a}^{(ij)}|\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}})} }

We can now relate these reconstructed entries to the data we have (likelihood function). Firstly, the vector \boldsymbol{a}^{(ij)} as mentioned has some associated input data. We have some experimental data for ij which is denoted by \boldsymbol{y}^{(ij)} which has input data \boldsymbol{x}\_{y,T}^{(ij)}, \boldsymbol{x}\_{y,M}^{(ij)}. That is we can look at \boldsymbol{a}^{(ij)} as being interpolated data from the experimental data \boldsymbol{y}^{(ij)} which we use to learn the lower rank matrices \boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}} but in a single framework. The likelihood can thus be described by the predictive equations of a GP, P(\boldsymbol{y}^{(ij)}|\boldsymbol{a}^{(ij)}). Or across all the known dataset:

P(\boldsymbol{y}|\boldsymbol{\mathcal{A}}) = \prod\_{i=1}^{N} {\prod\_{j=i+1}^{N} { {(P(\boldsymbol{y}^{(ij)}|\boldsymbol{a}^{(ij)}))^{I\_{ij}}} }}

Where: I\_{ij} is an indicator function equals to 1 if dataset ij is known, zero otherwise. The full bayesian model is thus:

P(\boldsymbol{\mathcal{A}}, \boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}}|\boldsymbol{y}) \propto P(\boldsymbol{y}|\boldsymbol{\mathcal{A}}) P(\boldsymbol{\mathcal{A}}|\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}}) P(\boldsymbol{\mathcal{U}}) P(\boldsymbol{\mathcal{V}})

I have the following stan for this implementation for this code

```no-highlight
functions {
    // kernel for composition
    matrix Kx(vector x1, vector x2, int order) {
        int N = rows(x1);
        int M = rows(x2);
        matrix[N, order+1] X1;
        matrix[M, order+1] X2;
        for (i in 1:order) {
            X1[:,i] = x1 .^(order+2-i) - x1;
            X2[:,i] = x2 .^(order+2-i) - x2;
        }
        X1[:,order+1] = 1e-1 * x1 .* sqrt(1-x1) .* exp(x1);
        X2[:,order+1] = 1e-1 * x2 .* sqrt(1-x2) .* exp(x2);
        return X1 * X2';
    }

    // kernel for Temperature
    matrix KT(vector T1, vector T2) {
        int N = rows(T1);
        int M = rows(T2);
        matrix[N, 4] TT1 = append_col(append_col(append_col(rep_vector(1.0, N), T1), 1e-2* T1.^2), 1e-4* T1.^3);
        matrix[M, 4] TT2 = append_col(append_col(append_col(rep_vector(1.0, M), T2), 1e-2* T2.^2), 1e-4* T2.^3);

        return TT1 * TT2';
    }

    // Combined kernel
    matrix K(vector x1, vector x2, vector T1, vector T2, int order) {
        return Kx(x1, x2, order) .* KT(T1, T2); 
    }

    // Reduce sum for parallel computing of log densities
    // Scaled Feature matrices priors
    real ps_F_raw(array[] int M_slice, int start, int end, array[] matrix F_raw) {
        real F_raw_target = 0;
        int M = num_elements(M_slice); // number of a_MC in the parallel space
        for (m in 1:M) {
            F_raw_target += std_normal_lpdf(to_vector(F_raw[2*M_slice[m]-1,:,:])); // Prior for scaled U matrices
            F_raw_target += std_normal_lpdf(to_vector(F_raw[2*M_slice[m],:,:])); // Prior for scaled V matrices
        }
        return F_raw_target;
    }

    // Likelihood and smoothed excess enthalpy entries
    real ps_a_MC_y(array[] int N_slice, int start, int end, matrix y_MC, array[,] matrix F_raw, vector v_ARD, 
        int N_known, int N_unknown, int N_MC, int N_T, int M, array[,] int Idx_all, 
        array[] int N_points, vector y1, matrix K_y, vector v_D, matrix K_MC, vector x1, 
        vector T1, vector x2, vector T2, real error, int order, matrix L_MC) {
            int N = num_elements(N_slice); // number of a_MC in the parallel space
            matrix[N, N_MC] a_MC;
            real all_target = 0;

            for (i in 1:N) {
                int counter = 1;
                for (t in 1:N_T) {
                    // composition up until just before 0.5
                    for (m in 1:M-1) {
                        a_MC[i, counter] = F_raw[t,m*2-1,:,Idx_all[N_slice[i],1]]' * diag_matrix(v_ARD) * F_raw[t,m*2,:,Idx_all[N_slice[i],2]];
                        counter += 1;
                    }
                    // composition of 0.5
                    a_MC[i, counter] = F_raw[t,2*M-1,:,Idx_all[N_slice[i],1]]' * diag_matrix(v_ARD) * F_raw[t,2*M-1,:,Idx_all[N_slice[i],2]];
                    counter += 1;
                    // composition from just after 0.5
                    for (m in 1:M-1) {
                        a_MC[i, counter] = F_raw[t,(M-m)*2-1,:,Idx_all[N_slice[i],2]]' * diag_matrix(v_ARD) * F_raw[t,(M-m)*2,:,Idx_all[N_slice[i],1]];
                        counter += 1;
                    }
                }

                all_target += multi_normal_cholesky_lpdf( y_MC[N_slice[i], :]' | a_MC[i,:]', L_MC);

                if (N_slice[i]<=N_known) {
                    matrix[N_points[N_slice[i]], N_MC] K_y_MC = K(x1[sum(N_points[:N_slice[i]-1])+1:sum(N_points[:N_slice[i]])], x2, T1[sum(N_points[:N_slice[i]-1])+1:sum(N_points[:N_slice[i]])], T2, order); // GP kernel for the cross covaraince
                    matrix[N_points[N_slice[i]], N_MC] stable_cov_y = mdivide_right(K_y_MC, (cholesky_decompose(K_MC))') ; // stable decomposition the K_(y-MC)*K_(MC)^-1*K_(y-MC)^T
                    matrix[N_points[N_slice[i]], N_points[N_slice[i]]] L_y = cholesky_decompose(K_y[sum(N_points[:N_slice[i]-1])+1:sum(N_points[:N_slice[i]]), sum(N_points[:N_slice[i]-1])+1:sum(N_points[:N_slice[i]])] + diag_matrix(error*abs(y1[sum(N_points[:N_slice[i]-1])+1:sum(N_points[:N_slice[i]])])+v_D[N_slice[i]]) - stable_cov_y*stable_cov_y'); //covariance for the likelihood function
                    vector[N_points[N_slice[i]]] mean_y = mdivide_right_spd(K_y_MC, K_MC) * y_MC[N_slice[i], :]'; //mean for the likleihood function
                    all_target += multi_normal_cholesky_lpdf(y1[sum(N_points[:N_slice[i]-1])+1:sum(N_points[:N_slice[i]])] | mean_y, L_y);
                }
            }
            return all_target;
    }
    
}

data {
    int N_known; // number of known mixtures
    int N_unknown; // number of unknown mixtures
    array[N_known] int N_points; // number of experimental points per known dataset
    int order; // order of the compositional polynomial
    vector[sum(N_points)] x1; // experimental composition
    vector[sum(N_points)] T1; // experimental temperatures
    vector[sum(N_points)] y1; // experimental excess enthalpy
    int N_T; // number of interpolated temperatures
    int N_C; // number of interpolated compositions
    vector[N_T] T2_int; // unique temperatures to interpolate
    vector[N_C] x2_int; // unique compositions to interpolate

    int N; // number of components
    int D; // rank of feature matrices
    array[N_known, 2] int Idx_known; // indices (row, column) of known datasets
    array[N_unknown, 2] int Idx_unknown; // indices (row, column) of unknown datasets 
    vector<lower=0>[N_known+N_unknown] v_D; // data-model mistamatch error
    real<lower=0> v_MC; // error between smoothed y_MC values and reconstructed entries

    real<lower=0> jitter; // jitter to stabalize inverse calculations
}

transformed data {
    real error=0.01;
    int M = (N_C + 1) %/% 2; // interger division to get the number of U matrices
    int N_MC = N_C*N_T; // overall number of interpolated datapoints per dataset
    vector[N_MC] x2; // concatnated vector of x2_int
    vector[N_MC] T2; // concatenated vector of T2_int
    matrix[N_MC, N_MC] K_MC; // kernel for the interpolated data
    matrix[sum(N_points), sum(N_points)] K_y; // kernel for the experimental data
    array[N_known+N_unknown] int N_slice;
    array[M-1] int M_slice;
    matrix[N_MC, N_MC] K_MC_inv;
    matrix[N_MC, N_MC] cov_y_MC;
    matrix[N_MC, N_MC] L_MC; // cholesky decomposition of the covariance of y_MC

    for (i in 1:N_known+N_unknown) {
        N_slice[i] = i;
    }

    for (i in 1:M-1) {
        M_slice[i] = i;
    }

    // Assign MC input vectors for temperature and composition
    for (i in 1:N_T) {
        x2[(i-1)*N_C+1:i*N_C] = x2_int;
        T2[(i-1)*N_C+1:i*N_C] = rep_vector(T2_int[i],N_C);
    } 

    // Assign MC kernel
    K_MC = K(x2, x2, T2, T2, order) + diag_matrix(rep_vector(jitter, N_MC));

    // Assign experimental data kernel
    for (i in 1:N_known) {
        K_y[sum(N_points[:i-1])+1:sum(N_points[:i]), sum(N_points[:i-1])+1:sum(N_points[:i])] = K(x1[sum(N_points[:i-1])+1:sum(N_points[:i])], x1[sum(N_points[:i-1])+1:sum(N_points[:i])], T1[sum(N_points[:i-1])+1:sum(N_points[:i])], T1[sum(N_points[:i-1])+1:sum(N_points[:i])], order);
    }

    // Compute cholesky decomposition of the cov of y_MC
    K_MC_inv = inverse_spd(K_MC);
    K_MC_inv = (K_MC_inv+K_MC_inv')/2; //ensure symmetricness
    cov_y_MC = inverse_spd(K_MC_inv+ diag_matrix(rep_vector(v_MC^-1,N_MC)));
    L_MC = cholesky_decompose(cov_y_MC);
}

parameters {
    // MC parameters
    array[N_T, M*2-1] matrix[D,N] F_raw; // scaled feature matrices; M U matrices, and M-1 V matrices
    vector<lower=0>[D] v_ARD; // variance of the ARD framework
    real<lower=0> scale; // scale parameter dictating strenght of ARD effect 

    // Smoothed GP parameters
    matrix[N_known+N_unknown, N_MC] y_MC; // smoothed excess enthalpy values

    // Data parameters
    //vector<lower=0>[N_known+N_unknown] v_D; // variance for the data-model mismatch.
}

model {
    int grainsize_F_raw = 1;
    int grainsize_y_MC = 1;

    // prior on scale
    scale ~ gamma(2,1);

    // prior on v_ARD
    v_ARD ~ exponential(scale);

    // Prior on the data-model mismatch
    //v_D ~ exponential(1);

    // ARD priors for feature matrices
    for (t in 1:N_T) {
        target += reduce_sum(ps_F_raw, M_slice, grainsize_F_raw, F_raw[t,:2*M-1]);
        to_vector(F_raw[t,2*M-1,:,:]) ~ std_normal(); // Prior for scaled U matrix at x=0.5
    }

    target += reduce_sum(ps_a_MC_y, N_slice, grainsize_y_MC, y_MC, F_raw, v_ARD, 
        N_known, N_unknown, N_MC, N_T, M, append_array(Idx_known, Idx_unknown), 
        N_points, y1, K_y, v_D, K_MC, x1, 
        T1, x2, T2, error, order, L_MC);
}

```

1. I wanted to know if there is any way I can optimize this code further?

2. I have used the `reduce_sum` functions to specify the priors and likelihood on the model as in the code provided. I have tested this code and found that for the normal code (using just for loops and not `reduce_sum`) optimization took 30min - 1 hour, which has been reduced to about 14 min per 1000 optimization iterations. I have played around by changing the number of stan threads used from 4 up until 7 but this does not seem to have any influence on the run time (I have an 8 core machine but will run this code on a cluster with significantly more resources, but I do not want to use more cores if it will not help with faster runtimes). Does this occur due to the overhead costs being significantly more when we use more cores? I have also tried to change the grainsize of the 2 `reduce_sum` functions but these also did not seem to have a significant effect on the runtime. Any suggestions would be helpful.

3. Also with the normal code I ran (no `reduce_sum`) the code took about 2.5 GB of ram but with the new model it takes about 600 MB. This is something I did not expect to happen and the code is being ran in parallel thus it should’ve used more resources? But then again in my original code I specified a matrix in the model block `matrix[N_known+N_unknown, N_MC] a_MC;` which probably lead to the reduction in memory usage as I now specify this matrix in the `reduce_sum` function, as an `matrix[num_elements(N_slice), N_MC] a_MC`; which could lead to the reduction in memory usage?

For reference for the current dataset `int N_MC = 57`, `int N_known = 102`, `int N_unknown = 58`, `int N=17`, `int D=10`, `int N_T = 3`, `int N_C=19`, `N_points`(number of experimental datapoints per dataset ij) ranges from 7-100. These values will becomes larger when I consider more data.

Also I am using cmdstanpy 1.2.1 and cmdstan 2.34.0

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [March 18, 2024, 1:46pm UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492/2 "2024-03-18T13:46:24Z")

</div>

Hi,  
I only have partial answers as I don’t fully understand the model, but since nobody else answered this, I’ll give it a go.

> [@Garren\_Hermanus](#):
>
> I wanted to know if there is any way I can optimize this code further?

Just to state the possibly obvious: Stan performance is heavily influenced by the geometry of the posterior. Improving the geometry can thus lead to substantial speedups (because you reduce treedepth necessary to reach U-turn, so you compute less steps per iteration). So checking some pairs plots etc. to detect signs of possibly problematic geometry is likely a good time investment.

I also can’t imagine debugging such a huge model, so I would definitely try decomposing it into a set of smaller components that can be tested and optimized separately and combined them only once the components work good in isolation (e.g. with SBC, following the [Small model implementation workflow • SBC](https://hyunjimoon.github.io/SBC/articles/small_model_workflow.html))

One possible optimizationin would be to use Hilbert-space basis approximation for the Gaussian process - see [Gaussian process demonstration with Stan](https://avehtari.github.io/casestudies/Motorcycle/motorcycle.html) but note that you’ll have to derive the basis specifically for your kernel of interest.

There might be tiny bits of improvement in using compound opeartion functions like `add_diag` or `dot_self` instead of raw matrix operations whenever possible. I’d say that constructs like `X' * diag_matrix(v_ARD) * Y` should be equivalent to `X' * Y * v_ARD` - which would be a bit faster if I am right.

You are also repeating the computation of `cholesky_decompose(K_MC)` within a loop in `ps_a_MC_y`.

I _think_ (not 100% sure) that if you are doing `cholesky_decompose` yourself on the covariance matrix, you may as well just run `multi_normal` directly with the full matrix and this might be a bit faster (at least is not gonna save the decomposition on the autodiff stack).

> [@Garren\_Hermanus](#):
>
> I have used the `reduce_sum` functions to specify the priors and likelihood on the model as in the code provided. I have tested this code and found that for the normal code (using just for loops and not `reduce_sum`) optimization took 30min - 1 hour, which has been reduced to about 14 min per 1000 optimization iterations. I have played around by changing the number of stan threads used from 4 up until 7 but this does not seem to have any influence on the run time (I have an 8 core machine but will run this code on a cluster with significantly more resources, but I do not want to use more cores if it will not help with faster runtimes). Does this occur due to the overhead costs being significantly more when we use more cores? I have also tried to change the grainsize of the 2 `reduce_sum` functions but these also did not seem to have a significant effect on the runtime. Any suggestions would be helpful.

Within-chain paralellization tends to have substantial overhead in moving values between threads. So yes, increasing the number of threads available may quickly reach diminishing returns. In that sense your situation appears typical.

> [@Garren\_Hermanus](#):
>
> Also with the normal code I ran (no `reduce_sum`) the code took about 2.5 GB of ram but with the new model it takes about 600 MB. This is something I did not expect to happen and the code is being ran in parallel thus it should’ve used more resources?

This is somewhat expected - the difference is that with `reduce_sum`, the backward pass of the autodiff is computed separately for each slice and then only the final value and its gradients are passed back to the main thread. Without `reduce_sum`, Stan will build the autodiff stack for the whole model and then do a single reverse pass once the `model` block completes. For this reason, people sometimes report reduced memory usage (and even performance gains due to e.g. better memory locality) when using `reduce_sum` even when just a single core is used.

Hope that helps you move forward!

---

<div class="post-metadata">

**Author:** ![Garren\_Hermanus](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/garren_hermanus/32/17356_2.png) [@Garren\_Hermanus](https://discourse.mc-stan.org/u/Garren_Hermanus)\
**Post date:** [March 18, 2024, 2:37pm UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492/3 "2024-03-18T14:37:49Z")

</div>

Thanks for the response @martinmodrak it is much appreciated!

> [@martinmodrak](#):
>
> I also can’t imagine debugging such a huge model, so I would definitely try decomposing it into a set of smaller components that can be tested and optimized separately and combined them only once the components work good in isolation (e.g. with SBC, following the [Small model implementation workflow • SBC](https://hyunjimoon.github.io/SBC/articles/small_model_workflow.html))

I have tested the model out and it works only when we provide initial values for \boldsymbol{\mathcal{U}} and \boldsymbol{\mathcal{V}} such that \boldsymbol{a}^{(ij)}\approx (\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}})\_{(ij)} where (\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{V}})\_{(ij)} is the product of the the matrix factorization where we only extract entries across ij and ji. When these are not provided the model tends to learn \boldsymbol{\mathcal{U}} and \boldsymbol{\mathcal{V}} to be zeros such that \boldsymbol{a}^{(ij)} is only the interpolated data from the experimental values. This would indicate to me that the posterior is most likely multimodal, and given that the matrix factorization is unidentifiable contributes more towards modelling difficulties. But I will check this post and see if it might be useful for the application of my model.

> [@martinmodrak](#):
>
> There might be tiny bits of improvement in using compound opeartion functions like `add_diag` or `dot_self` instead of raw matrix operations whenever possible

I’ll have a look at the documentation for this again, thanks.

> [@martinmodrak](#):
>
> I’d say that constructs like `X' * diag_matrix(v_ARD) * Y` should be equivalent to `X' * Y * v_ARD` - which would be a bit faster if I am right.

I do not think this would be valid since X, Y are both D dimensional vectors, or if we consider the full matrices these are N\times D dimensional matrices and `diag_matrix(v_ARD)` is a D\times D diagonal matrix. That is the latent dimensions of Y and X is scaled by v\_ARD, similar to SVD. A similar construct would probably be to dot multiple and then dot product. That is `dot_product(X.*v_ARD, Y)`. I’ll test this out and see if it makes a difference or not.

> [@martinmodrak](#):
>
> You are also repeating the computation of `cholesky_decompose(K_MC)` within a loop in `ps_a_MC_y`.

Thanks, I made some adjustments to the code where previously I incorporated noise into the kernel matrix.

> [@martinmodrak](#):
>
> I _think_ (not 100% sure) that if you are doing `cholesky_decompose` yourself on the covariance matrix, you may as well just run `multi_normal` directly with the full matrix and this might be a bit faster (at least is not gonna save the decomposition on the autodiff stack).

I followed the user guide on this, which stated that using the `multi_normal_cholesky` is faster than `multi_normal`. This may have been an outdated version of the documentation which I read as I do not see this in the current documentation. I’ll see if this makes much of a difference.

> [@martinmodrak](#):
>
> Within-chain paralellization tends to have substantial overhead in moving values between threads. So yes, increasing the number of threads available may quickly reach diminishing returns. In that sense your situation appears typical.

I see, thanks. I’ll just stick to 4 cores for this

Thanks again for your response!

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [March 18, 2024, 8:42pm UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492/4 "2024-03-18T20:42:54Z")

</div>

> [@Garren\_Hermanus](#):
>
> I do not think this would be valid since XXX, YYY are both D dimensional vectors, or if we consider the full matrices these are N\times DN×DN\times D dimensional matrices and `diag_matrix(v_ARD)` is a D\times DD×DD\times D diagonal matrix. That is the latent dimensions of YYY and XXX is scaled by v\_ARD, similar to SVD. A similar construct would probably be to dot multiple and then dot product. That is `dot_product(X.*v_ARD, Y)`. I’ll test this out and see if it makes a difference or not.

You are correct, I did think that `v_ARD` is a scalar, but that was a misunderstanding.

> [@Garren\_Hermanus](#):
>
> This would indicate to me that the posterior is most likely multimodal, and given that the matrix factorization is unidentifiable contributes more towards modelling difficulties.

Yes, multimodality/non-identifiability can hurt your performance substantially.

Good luck going forward.

---

<div class="post-metadata">

**Author:** ![mitzimorris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mitzimorris/32/29_2.png) [@mitzimorris](https://discourse.mc-stan.org/u/mitzimorris)\
**Post date:** [March 19, 2024, 4:33pm UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492/5 "2024-03-19T16:33:43Z")

</div>

not an optimization, but since you’re using CmdStanPy, have you tried running Pathfinder first to improve the inits?

example here: [Variational Inference using Pathfinder — CmdStanPy 1.2.1 documentation](https://mc-stan.org/cmdstanpy/users-guide/examples/Pathfinder.html#Pathfinders-as-initialization-for-the-MCMC-sampler)

---

<div class="post-metadata">

**Author:** ![Garren\_Hermanus](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/garren_hermanus/32/17356_2.png) [@Garren\_Hermanus](https://discourse.mc-stan.org/u/Garren_Hermanus)\
**Post date:** [March 19, 2024, 4:46pm UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492/6 "2024-03-19T16:46:45Z")

</div>

I did not know something like this existed. I’ll have a look at this. I have tried a bunch of different initialization such as running MAP and VI from random points and using the mean from VI or optimized params from MAP. I do not know if this is a similar approach to what pathfinder is?

---

<div class="post-metadata">

**Author:** ![Garren\_Hermanus](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/garren_hermanus/32/17356_2.png) [@Garren\_Hermanus](https://discourse.mc-stan.org/u/Garren_Hermanus)\
**Post date:** [May 18, 2024, 3:07am UTC](https://discourse.mc-stan.org/t/optimizing-stan-code-runtime-for-matrix-factorization/34492/7 "2024-05-18T03:07:06Z")

</div>

Coming back to this, I have made an intermediate reparameterization.  
Changing the following multivariate normal to a standard normal using the cholesky decomposition.

> [@Garren\_Hermanus](#):
>
> `all_target += multi_normal_cholesky_lpdf( y_MC[N_slice[i], :]' | a_MC[i,:]', L_MC);`

The speed up times I am experiencing is really significant. This initially took ±2 hours to draw 200 samples (100 warmup + 100 samples). After the reparameterization it takes ±1 hour to draw 2000 samples (1000 warmup + 1000 samples). I am not specifying the original parameter in the transformed parameters block, but rather using the raw parameter in the `partial_sum` likelihood function.

Q: I just want to know why this speedup is so significant when moving to a non-cenetered case?