fitlm warns of singular design matrix when run in parfor loop, but not in a serial for loop (same data)

I wished to estimate the distributions of some linear regression beta values using a bootstrap resampling. To save time, I ran the fitlm beta estimation of bootstrap samples in a parfor loop. To do so, I did the following:

1. Set random seed for reproducibility.
2. Generated bootstrap replicate indices (indices of original data sampled with replacement)
3. Passed original data table and a bootstrap replicate's indices to a worker using parfor.
4. Constructed the bootstrapped data on the worker using the original data table and bootstrap indices
5. Ran fitlm on the worker and saved the beta values from the linear regression

When I run this code, I get the following warning: "Warning: Regression design matrix is rank deficient to within machine precision." This warning appears for many (but not all) of my bootstrap replicates.

However, when I run the above process serially, I do not get the rack-deficiency warning. The bootstrapped data should be the same between serial and parallel cases as I am setting the random seed before generating the boostrap replicates; I simply switch the parfor to for. I am not suppressing warnings in either case.

For what it is worth, my data table contains both continuous and categorical variables. When I limit my model to continuous variables, I do not get the rank-deficiency warning when I run the code in parallel. This warning is only generated (and only in the parallel case) when I include categorical variables in my model. I have roughly 13 million rows in my data, so unfortunately I can't share it here for error replication purposes (>1GB in size). However, with this many data points it seems unlikely a bootstrap replicate would generate a rank deficient design matrix by chance - let alone many rank deficient matrices.

When I compare the beta values generated by parallel and serial fitlm execution, they do differ but only at the level of machine precision (mean absolute difference: 8.8e-16). So despite throwing the rank deficiency error in parallel, fitlm seems to be arriving at very similar beta estimates regardless of whether it's run in parallel or serially.

Any feedback on this would be greatly appreciated! I have not previously seen a case of parallel-specific warnings when using the parallel toolbox. I am speculating that the data table might be modified (compressed?) when it is passed to the workers, leading to problems with the design matrix. But this is just speculation. I didn't see anything in the documentation to reflect this.

Thanks!

Pseudo-code below:

rng('default');     % Set random seed
rng(1);

[~,bootsamp]=bootstrp(10000,[],1:height(bigDataTable)); % Generate 10k bootstrap replicates

betas_par = nan(10000,5);   % Variables to hold beta estimates
betas_ser = nan(10000,5);

% This code leads to rank-deficient design matrix warnings: 
parfor i = 1:10000
    bootDataTable = bigDataTable(bootsamp(:,i),:);  % Create boostrapped data on worker
    mdl = fitlm(bootDataTable, 'y~1+x1+x2+x3+x4');  % Fit betas on worker (x1-x4 are continuous or categorical vars)
    betas_par(i,:) = mdl.Coefficients.Estimate';    % Save beta coefficients
end

% This almost-identical code does not produce warnings despite running the same models with the same data:
for i = 1:10000
    bootDataTable = bigDataTable(bootsamp(:,i),:);  % Create boostrapped data
    mdl = fitlm(bootDataTable, 'y~1+x1+x2+x3+x4');  % Fit betas (x1-x4 are continuous or categorical vars)
    betas_ser(i,:) = mdl.Coefficients.Estimate';    % Save beta coefficients
end

[Bruno Luong]

"The bootstrapped data should be the same between serial and parallel cases as I am setting the random seed before generating the boostrap replicates; I simply switch the parfor to for. I am not suppressing warnings in either case."

Not sure if it's matter but the seed setting before the parfor loop has no effect.

Maybe your random data have defect because of this.

[James Hengenius]

Hi Bruno - Thanks for the feedback! While I see what you mean that setting the random seed shouldn't affect RNG calls on the workers, I call the RNG outside of the parfor loop (using the bootstrp function) and then pass the results to the parfor workers. So I don't think this is what's causing the warning. (Though I am not 100% sure.)

[Walter Roberson]

I wonder if it uses random numbers for anything? Even just initialization of search position?

[James Hengenius]

Hi Walter - do you mean fitlm? Possibly the search is initialized differently, though I think before fitlm begins the search it checks the design matrix using classreg.regr.CompactTermsRegression/checkDesignRank. This checkDesignRank is throwing the warning in the parallel case but not the serial case - before the search begins as far as I can tell.

[Bruno Luong]

May be one of the function you use set the seed. It's very dangerous when runing in parallel.

You might check (in regular for) who calls rng()

[James Hengenius]

This documentation confirms that different random number generators are used on workers vs. client: https://www.mathworks.com/help/parallel-computing/control-random-number-streams-on-workers.html

If I'm reading this correctly, it appears setting a seed on the client should not effect the workers. However, as Walter pointed out, this difference in random generators might be affecting the execution of fitlm on the workers.

[Bruno Luong]

Especially if rng('shuffle') is called, then the random on dofferent workers are the same, which is not good.

Again I recommend you should chase who calls rng.

" it appears setting a seed on the client should not effect the workers."

This is sure as I told you earlier.

[Walter Roberson]

parallel workers use different streams of the random generation.

The question now is whether fitlm needs random numbers for anything.

[Bruno Luong]

The bootstrap data are constructed using random selection, as I understand.

[James Hengenius]

I tried setting the workers' RNG to Twister (the client's default RNG) to see if the choice of RNG on the workers was affecting the output:

% This code leads to rank-deficient design matrix warnings: 
parfor i = 1:10000
    
    rng(0,'Twister') % Set the RNG to client default
    
    bootDataTable = bigDataTable(bootsamp(:,i),:);  % Create boostrapped data on worker
    mdl = fitlm(bootDataTable, 'y~1+x1+x2+x3+x4');  % Fit betas on worker (x1-x4 are continuous or categorical vars)
    betas_par(i,:) = mdl.Coefficients.Estimate';    % Save beta coefficients
end

However, the same rank-deficiency warning occurred regardless of whether the workers' RNG was Twister (client default) or Threefry (worker default).

[Bruno Luong]

TMW recommends using different randstream for different worker https://uk.mathworks.com/help/matlab/ref/randstream.html

[Edric Ellis]

One other thing to try - call

maxNumCompThreads(1)

in your client MATLAB before running the for loop version. This causes it to run the same way that the workers will run, and I have a recollection that some operations can be sensitive to that.

[Walter Roberson]

It looks like fitlm might not use random numbers.

load carsmall
tbl = table(Weight,Acceleration,MPG,'VariableNames',{'Weight','Acceleration','MPG'});
rng(12345)

start_rng = rng()

start_rng = struct with fields:
     Type: 'twister'
     Seed: 12345
    State: [625×1 uint32]

lm = fitlm(tbl,'MPG~Weight+Acceleration')

lm = 
Linear regression model:
    MPG ~ 1 + Weight + Acceleration

Estimated Coefficients:
                     Estimate         SE         tStat       pValue  
                    __________    __________    _______    __________

    (Intercept)         45.155        3.4659     13.028    1.6266e-22
    Weight          -0.0082475    0.00059836    -13.783    5.3165e-24
    Acceleration       0.19694       0.14743     1.3359       0.18493


Number of observations: 94, Error degrees of freedom: 91
Root Mean Squared Error: 4.12
R-squared: 0.743,  Adjusted R-Squared: 0.738
F-statistic vs. constant model: 132, p-value = 1.38e-27

end_rng = rng();

isequal(start_rng, end_rng)

ans = logical
   1

[Walter Roberson]

as an experiment, try asking to fit the same data on multiple workers. If you can track down which i value is giving the warning in your existing code, try that value on multiple workers

[James Hengenius]

So I think I found the reason why this warning is occurring. But I don't understand why the warning is only occurring on workers in my original code.

One of the categorical variables I am including is race. One of the possible values for this variable is Other. (Only one person in my sample has the value of Other.) So when I do bootstrap resampling, it is possible that a resulting bootstrap replicate will contain no rows with race == Other.

As a test, I excluded all rows with race == Other

tbl; % Is my full dataset
tbl_exclude = tbl; tbl_exclude(tbl_exclude.race=='Other')=[]; % Remove all rows with race==Other

Then I checked that the rows were in fact removed.

find(tbl_exclude.race=='Other')
ans =
 0×1 empty double column vector

Finally, I ran a test model in fitlm.

 fitlm(tbl_exclude, 'y~1+race')
Warning: Regression design matrix is rank deficient to within machine precision. 
> In classreg.regr/CompactTermsRegression/checkDesignRank (line 35)
In LinearModel.fit (line 1040)
In fitlm (line 121) 

ans = 

Linear regression model:
    y ~ 1 + race

Estimated Coefficients:
                                    Estimate       SE      tStat     pValue
                                   __________    ______    ______    ______

    (Intercept)                     6.521e+06    6781.8    961.55       0  
    race_Black/African American    4.0891e+05    6959.8    58.753       0  
    race_Other                              0         0       NaN     NaN  
    race_White                     3.7114e+05    6947.6    53.419       0  


Number of observations: 13801200, Error degrees of freedom: 13801197
Root Mean Squared Error: 3.98e+06
R-squared: 0.000255,  Adjusted R-Squared: 0.000255
F-statistic vs. constant model: 1.76e+03, p-value = 0

As you can see, it prduces the warning I was seeing. And even though the table I pass to fitlm doesn't contain any race==Other data, fitlm still tries to fit the beta for the variable race_Other.

I am guessing that the table tbl_exclude contains metadata from the original table tbl. This metadata might record that race==Other is a possibility even though it is no longer present in tbl_exclude. I looked in tbl_exclude.Properties for anything like this but didn't see it.

I am also unsure why my original code suppressed this warning when I ran the bootstrap serially.

EDIT TO ADD:

I found the function removecats (https://www.mathworks.com/help/matlab/ref/categorical.removecats.html) which removes unused categories from categorical arrays. After running

tbl_exclude.race = removecats(tbl_exclude.race)

I no longer get the rank deficient design matrix warning from fitlm.

This doesn't explain why the warning was suppressed in serial but not parallel execution, but it does identify the source of the warning. My bootstrap resample had a "ghost" category that was not represented in the data but was part of the categorical variables metadata. removecats cleaned this up.

Answer 1

[Bruno Luong]

Can you please try this code

[~,bootsamp]=bootstrp(10000,[],1:height(bigDataTable)); % Generate 10k bootstrap replicates
betas_par = nan(10000,5);   % Variables to hold beta estimates
EnvenBiggerDataTab = bigDataTable(bootsamp,:);
EnvenBiggerDataTab = reshape(EnvenBiggerDataTab, [size(bootsamp), size(bigDataTable,2)]);
EnvenBiggerDataTab = permure(EnvenBiggerDataTab,[1 3 2]);
% This code leads to rank-deficient design matrix warnings: 
parfor i = 1:10000
    bootDataTable = EnvenBiggerDataTab(:,:,i);  % Create boostrapped data on worker
    mdl = fitlm(bootDataTable, 'y~1+x1+x2+x3+x4');  % Fit betas on worker (x1-x4 are continuous or categorical vars)
    betas_par(i,:) = mdl.Coefficients.Estimate';    % Save beta coefficients
end

[James Hengenius]

Hi Bruno, I think I see what you're going for here. Unfortunately I don't have enough memory on my PC to make EvenBiggerDataTab. The original table is over a GB and increasing its size by 10,000 times is not feasible. I will try with a smaller bootstrap replicate number (maybe 3x) when I have some free time today.

My intuition based on your previous comments is that you want to avoid random number generator calls on the workers. However, my original code does not call the random number generator on the workers (as best I can tell based on Walter's comments). The variable bootsamp is generated on the client using the random number generator, but passing bootsamp to the workers should not require addiitional random number generator calls by the workers.

Is this what you are trying to test with this solution? Or am I missing something? Thanks!

[Bruno Luong]

I believe when you do

bootDataTable = bigDataTable(bootsamp(:,i),:)

inside parfor, there must be some complex exchanges between workers (I'm even surprise such syntax is accepted by MATLAB), or some concurrent access to a shared variable is going on. It is not clear to me how it handles internally.

I just suspect the problem could come from there.

May be you could try with a toy example rather than your final bootstrap to confirm this is an issue.

[Edric Ellis]

There is no worker-to-worker communication at all in parfor. In the statement:

bootDataTable = bigDataTable(bootsamp(:,i),:)

The variable bootsamp is a parfor sliced variable - only the relevant pieces of bootsamp are sent to each worker. The variable bigDataTable is a broadcast variable. The whole value is copied to each worker.

Answer 2

[Bruno Luong]

Could you try this:

C_bigDataTable = parallel.pool.Constant(bigDataTable);
% This code leads to rank-deficient design matrix warnings: 
parfor i = 1:size(bootsamp,2)
    bootDataTable = C_bigDataTable.Value(bootsamp(:,i),:);  % Create boostrapped data on worker
    mdl = fitlm(bootDataTable, 'y~1+x1+x2+x3+x4');  % Fit betas on worker (x1-x4 are continuous or categorical vars)
    betas_par(i,:) = mdl.Coefficients.Estimate';    % Save beta coefficients
end

[James Hengenius]

Hi Bruno, I gave this a try and made my bigDataTable a constant on the workers. Unfortunately, this had no effect. The parallel execution of fitlm still produces a rank deficient design matrix warning even with this setup.

[Bruno Luong]

I run out of idea.

If you can reproduce a small example may be you should report officially as a potential bug.