Coverage for src / sdynpy / signal_processing / sdynpy_lrm.py: 10%
228 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1import numpy as np
2import scipy as sp
3import math as mt
4try:
5 from joblib import Parallel, delayed
6except ImportError:
7 Parallel = None
8 delayed = None
9try:
10 import psutil
11except ImportError:
12 psutil = None
13import warnings
14from time import time
16def frf_local_model(references,responses,abscissa,f_out=None,bandwidth=None,
17 transient=True,modelset=[[1,1,0],[1,1,1],[1,1,2]],
18 export_ratio=.1,max_parallel=512,print_time=True):
19 """
20 Local modeling FRF estimator for MIMO systems.
21 This function can be called directly or using
22 frf.timedata2frf(...,method='LRM') or frf.fft2frf(...,method='LRM').
24 Parameters
25 ----------
26 references : ndarray (N_in, N_freqs) or (N_in, 1, N_freqs)
27 REFERENCES FREQUENCY-DOMAIN DATA.
28 responses : ndarray, shape (N_out, N_freqs) or (N_out, 1, N_freqs)
29 RESPONSES FREQUENCY-DOMAIN DATA.
30 abcissa : ndarray (N_freqs,)
31 INPUT FREQUENCY VECTOR.
33 f_out: ndarray (N_freqs_out,), optional
34 Output frequency vector. Finer resolution increases computational cost,
35 but only to a finite limit. Compuatational cost will not increase beyond
36 f_out = f_in, but large interpolation may result in a large export_ratio.
37 The default is f_in
38 bandwidth: float, optional
39 Local model estimation bandwidth in units of f_in. Larger values yield
40 better noise filtering but risk underfitting and take longer to run.
41 transient : boolean, optional
42 whetheter to include transient estimation in local models. Recommend False
43 for impact testing or burst random, True otherwise. The default is True.
44 modelset : iterable of numpy (3,) arrays, optional
45 Arrays contain [num. order, trans. order, denom. order].
46 Recommend increasing numerator order for nonlinear structures
47 and complicated FRFs. The default is [[1,1,0],[1,1,1],[1,1,2]].
48 export_ratio : f;pat, optional
49 (centered) proportion of local bands that can be stored and exported.
50 Larger values risk capturing volatile end behavior but yield faster
51 computation times. The default is 0.1.
52 max_parallel: int, optional
53 for very large systems with high frequency resolution and large bandwidths,
54 system RAM can be a constraint. The default number of parallel processes
55 is the number of CPU threads. If this causes freezing, max_parallel can
56 be set to reduce the allowed number of parallel threads. Also, smaller
57 pools may initialize faster. The default is 512.
58 print_time: boolean, optional
59 print remaining time if more than 10 seconds
61 Returns
62 -------
63 f_out: ndarray (N_freqs_out,)
64 estimated FRM frequency vector
65 H: ndarray (N_freqs_out, N_out, N_in)
66 frequency response matrix
67 model_data: dict
68 candidate model info and selected model at each frequency in f_out
71 Notes
72 -------
73 To improve noise and transient removal, recommend increasing bandwidth
75 To reduce computation time, recommend decreasing f_out resolution,
76 increasing export_ratio, or decreasing bandwidth. For very large systems,
77 if freezing occurs, reduce max_parallel.
79 For more information, see "A practitioner's guide to local FRF estimation",
80 K. Coletti. This function uses the MISO parameterization
82 """
84 """
85 NOTES FOR MAINTAINER
86 Outline of this function
87 -------
88 (1) f_in and f_out are parsed to generate an array of bands at which to
89 estimate local models (bands) and a list of query points in each
90 model at which to export solutions (query)
91 (2) based on the model candidate set and problem size, functional
92 parameterizations are generated which return linearized coefficient
93 matrices (for solving) and the FRF matrix (for export)
94 (3) computation time is estimated, and the LSQ algorithm is selected.
95 Solving model parameters using interative LLSQ is the majority of
96 computation time.
97 (4) in parallel, local models are solved for each model in modelset and
98 each frequency band. FRFs are exported. CTRL-F 'results = Parallel'
99 (5) the winning model and corresponding FRM is returned in each band
101 Notes on variable names
102 -------
103 U and Y are input and output data
104 A,B,D are FRF numerator, transient numerator, and denom., respectively
105 H is the frequency response matrix.
106 LR is a function that produces linear coefficients (L) and remainder (R) for the local models
107 See "A Practitioner’s Guide to Local FRF Estimation", K. Coletti for more information
108 """
110 if psutil is None and Parallel is None:
111 raise ValueError('frf_local_model requires modules psutil and joblib.\nPlease install these modules prior to using this function.')
112 elif psutil is None:
113 raise ValueError('frf_local_model requires the module psutil.\nPlease install this module prior to using this function.')
114 elif Parallel is None:
115 raise ValueError('frf_local_model requires the module joblib.\nPlease install this module prior to using this function.')
117 #%% Process arguments
118 U = references
119 Y = responses
120 f_in = abscissa
121 del references; del responses; del abscissa
123 df = float(f_in[1]-f_in[0]) # frequency step size
124 if bandwidth is None: # local model bandwidth
125 bandwidth = max( [min([300*df,40]) , max([30*df,10])] )
126 # bounded below by the larger of 30 points and 10 Hz,
127 # above by the smaller of 300 points and 40 Hz
128 if f_out is None: # spacing between center frequencies for estimation
129 f_out = f_in
131 # Default parameters, hidden from user
132 NitSK = 50 # max SK iterations
133 conv_tol = .005 # SK parameter convergence tolerance
135 #%% Calculate general helper vars
136 Nfreq = f_in.size # number of frequencies
137 Nin = U.shape[0]
138 Nout = Y.shape[0]
139 Nmod = len(modelset)
140 U = U.reshape(Nin,Nfreq,order='F') # no realizations/frames needed
141 Y = Y.reshape(Nout,Nfreq,order='F')
142 f_in = f_in.reshape(-1,order='F')
143 f_out = f_out.reshape(-1,order='F')
145 if not Y.shape[-1] == U.shape[-1] == Nfreq:
146 raise Exception('input arrays are not correctly sized')
148 #%% Get estimation bands and query points (STEP 1)
149 # Estimation bands contain the frequency bins used for parameter estimation,
150 # and query points are the frequencies at which results are exported
151 Nband = __round_odd(bandwidth/df) # number of data points in band
152 N_half_band = int(Nband/2) # number of data points on either side of center
153 max_from_center = export_ratio/2*bandwidth # maximum query distance from center point, except at ends
155 bands,query = __parse_bands(N_half_band,max_from_center,f_in,f_out) # bands gives indices of f_in, and query gives r values for query. r=Hz/df.
156 Ninds = bands.shape[0]
158 #%% Initialize parameterizations (STEP 2)
159 LR,D,lsq_err,H = [[None]*Nmod,[None]*Nmod,[None]*Nmod,[None]*Nmod]
160 Npar = np.zeros(Nmod)
162 for k in range(Nmod-1,-1,-1):
163 if transient:
164 LR[k],D[k],lsq_err[k],H[k],Npar[k] = __parameterize(Nin,Nout,Nband,modelset[k])
165 else:
166 LR[k],D[k],lsq_err[k],H[k],Npar[k] = __parameterize_notrans(Nin,Nout,Nband,modelset[k])
168 if Npar[k]>=Nband*Nout: # remove infeasible models
169 warnings.warn('bandwidth too small for ' + str(modelset[k]) + ' model order')
170 modelset.pop(k)
171 Nmod -= 1
172 LR.pop(k);D.pop(k);lsq_err.pop(k);H.pop(k);Npar = np.delete(Npar,k)
174 dof = Nband*Nout-Npar
176 #%% Compute (STEP 3/4)
177 # Parallel function. # joblib can't parse lists (only arrays), so query is a parameter
178 def __solve_ind(k,inq,lstsq):
179 # Initialize. inq are radii of H to output in this band
180 H_multi_model = np.empty((Nout,Nin,Nmod,inq.size),dtype=complex) # axis (0) output (1) input (2) model (3) frequency
181 mdl = np.empty((Nmod,inq.size),dtype=complex)
182 Uk = U[:,bands[k]]
183 Yk = Y[:,bands[k]]
185 for kk in range(Nmod): # model loop
186 H_multi_model[:,:,kk,:],pstar = __solveLRM(Uk,Yk,LR[kk],D[kk],lsq_err[kk],H[kk],inq,conv_tol,NitSK,lstsq)
188 # Model selection - MDL (see Pintelon, IEEE, 2021)
189 mdl[kk,:] = pstar/dof[kk] * mt.exp( mt.log(2*Nband)*Npar[kk]/(dof[kk]-2) )
191 return H_multi_model,mdl
193 # Estimate time, choose algorithm - SEE APPENDIX NOTES 1 AND 2
194 memory = psutil.virtual_memory().total/1e9
195 n_parallel = min(psutil.cpu_count(),61,max_parallel)
196 if memory/n_parallel < 2 and Nin*Nout > 500:
197 warnings.warn('Many parallel processes relative to system RAM. For large problems, freezing may occur.')
198 def test_lstsq_time(alg):
199 start = time()
200 __solve_ind( round(Ninds/2),query[round(Ninds/2)],alg )
201 return time()-start
202 if Nin*Nout < 500:
203 time_numpy = test_lstsq_time(__lstsq_numpy)
204 else:
205 time_numpy = np.inf # numpy is very slow for large problems
206 time_scipy = test_lstsq_time(__lstsq_scipy)
208 if time_numpy < time_scipy:
209 lstsq = __lstsq_numpy
210 else:
211 lstsq = __lstsq_scipy
213 time_est = min(time_numpy,time_scipy)*Ninds/n_parallel
214 if time_est > 10 and print_time: # parallel initialization takes 5 seconds
215 print( 'Estimated time remaining: parallel init. + ' + str( round(1.1*time_est,1)) + ' seconds' )
217 # Compute and extract results into arrays
218 results = Parallel(n_jobs=n_parallel)(delayed(__solve_ind)(k,query[k],lstsq) for k in range(Ninds))
219 H_multi_model = np.concatenate([x[0] for x in results],axis=3) # axis (0) Nout (1) Nin (2) model (3) Nind
220 mdl = np.concatenate([x[1] for x in results],axis=1) # axis (0) model (1) Nind
222 #%% Select winning models and export (STEP 5)
223 msel = np.argmax(mdl,0) # axis. (0) freq. ind
224 H = H_multi_model[:,:,msel,np.arange(f_out.size)] # axis. (0) Nout (1) Nin (2) freq. ind
225 model_data = {'info': 'Model format is [H num. order, trans. num. order, denom. order].',
226 'modelset':modelset,
227 'model_selected':msel}
229 # Return results
230 return f_out, np.moveaxis(H,-1,0), model_data
232#%% Local model solver
233def __solveLRM(U,Y,LR,D,lsq_err,H,query,conv_tol,NitSK,lstsq):
234 """
235 Local model solver. This function occupies most of the FRF local modeling estimation
236 time. Inputs are functions that accept parameter values and output local model component
237 matrices used for least-squares solving and model exporting. Uses Sanathanan-Koerner
238 iterations to solve the nonlinear rational LSTSQ problem iteratively
239 """
240 # SK iteration preparations
241 L,R = LR(U,Y) # linearized transfer matrix and residual vector. Lin objective = D^{-1}*(L*t-R)
242 Dp = np.ones(Y.shape) # initial scaling matrix. Axis (0) Nout (1) Nband
243 t0s,lsq_min = [float('inf')]*2
245 # Snathanan-Koererner Iterations
246 for k in range(NitSK+5):
247 # Define adjusted LLSQ/GLSQ problem
248 X = L/Dp[:,:,np.newaxis] # left matrix in LLSQ stacked over r. Axis (0) Nout (1) Nband (2) Npar
249 q = R/Dp # residual vector stacked over r. Axis (0) Nout (1) Nband
251 # Solve and advance. Vast majority of compution time in all of get_frf
252 t0 = lstsq(X,q) # axis. (0) Nout (1) Npar
254 # Convergence
255 if abs((t0-t0s)/t0).mean() < conv_tol:
256 break # exit if convergence is reached
257 else:
258 Dp = D(t0) # store new adjustment matrix
259 t0s = t0
261 # Capture nonconvergence (in my experience, this is very rare)
262 if k > NitSK:
263 lsq_t0 = lsq_err(U,Y,t0)
264 if lsq_t0 < lsq_min:
265 lsq_min = lsq_t0
266 t1 = t0
268 t0=t1 if 't1' in locals() else t0
270 return H(t0,query),-lsq_err(U,Y,t0) # Hstar, pstar (log density)
272#%% Parameterizations. See "A Practitioner’s Guide to Local FRF Estimation", K. Coletti
273def __parameterize(Nin,Nout,Nband,order):
274 """
275 Generate parameterization functions for estimation including transient removal
276 Parameterizations are passed into __solveLRM to generate local models
277 """
278 # Components of f, because f can be iterated
279 # Radius and exponent matrices. r is radius, s power, and A/B/D component matrices
280 N_half_band = (Nband-1)/2
281 r = np.arange(-N_half_band,N_half_band+1)
282 rsA = r**np.arange(order[0]+1)[:,np.newaxis] # axis. (0) power (1) radius
283 rsB = r**np.arange(order[1]+1)[:,np.newaxis]
284 rsD = r**np.arange(order[2]+1)[:,np.newaxis]
286 # Define theta break points
287 bA = np.arange(Nout*Nin*(order[0]+1)) # transfer numerator
288 bB = np.arange(bA[-1]+1,bA[-1]+Nout*(order[1]+1)+1) # transient numerator
289 bD = np.arange(bB[-1]+1,bB[-1]+Nout*order[2]+1) # denominator (wrap Nout first, then s)
290 Npar = bB[-1]+Nout*order[2]+1
292 def lsq_err(U,Y,t):
293 t = t.reshape(-1,order='F') # transform from SO format to MO
295 # Build matrices
296 A = np.einsum('ij,jk->ik',t[bA].reshape(-1,order[0]+1,order='F'),rsA).reshape(Nout,Nin,-1,order='F')
297 B = np.einsum('ij,jk->ik',t[bB].reshape(-1,order[1]+1,order='F'),rsB) # axis (0) Nout (1) Nfreq
298 catmat = np.concatenate( (np.ones((Nout,1)),t[bD][:,np.newaxis]),0 )
299 D = np.einsum('ij,jk->ik',catmat.reshape(Nout,-1,order='F'),rsD) # axis. (1) Nout (2) Nfreq
301 # Compute sum of squared residuals
302 Ysim = (np.einsum('jkl,kl->jl',A,U) + B)/D
303 return (abs(Y-Ysim)**2).sum()
305 def H(t,r_inq):
306 t = t.reshape(-1,order='F') # transform from SO format to MO
308 # Radius and exponent matrices
309 rsA_inq = r_inq**np.arange(order[0]+1)[:,np.newaxis] # axis. (0) power (1) radius
310 rsD_inq = r_inq**np.arange(order[2]+1)[:,np.newaxis]
312 # Build matrices
313 A = np.einsum('ij,jk->ik',t[bA].reshape(-1,order[0]+1,order='F'),rsA_inq).reshape(Nout,Nin,-1,order='F') # axis. (1) Nout (2) Nin (3) Nfreq
314 temp = np.concatenate( (np.ones(Nout),t[bD]) )
315 D = np.einsum('ij,jk->ik',temp.reshape(Nout,-1,order='F'),rsD_inq) # axis. (1) Nout (2) Nfreq
317 return A/D[:,np.newaxis,:]
319 # Build linearized (SK) transfer matrix L(U,y)*t + R(U,Y) = Dp^(-1) (D(t)*Y - N(t)*U - M(t))
320 # Denominator matrix for SK iterations
321 def D(t):
322 return 1 + np.einsum( 'ij,lj->li',r[:,np.newaxis]**np.arange(1,order[2]+1) , t[:,np.arange(-order[2],0)] ) # axis (0) Nout (1) Nband
324 def LR(U,Y):
325 # Compute residual
326 R = -Y # axis (0) Nout (1) Nband
328 # Compute components
329 r1 = r[:,np.newaxis]
330 A1 = Y.T.reshape(-1,1,Nout,order='F')* (r1**np.arange(1,order[2]+1))[:,:,np.newaxis] # D(t)*Y. Axis. (1) r (2) s (3) Nout
331 A2 = (U.T[:,:,np.newaxis] * r1[:,:,np.newaxis] **np.arange(order[0]+1)).reshape(Nband,Nin*(order[0]+1),1,order='F') # A(t)*U. Axis. (1) r (2) Nin*s (3) Nout
332 A3 = ( r1**np.arange(order[1]+1) )[:,:,np.newaxis] # B(t). Axis. (1) r (2) s (3) empty
334 # Assemble into array
335 if order[-1]>0:
336 catmat = np.ones((1,1,Y.shape[0])) # A2 and A3 don't depend on Nout
337 L = np.concatenate((-A2*catmat,-A3*catmat,A1),1).transpose(2,0,1) # axis (0) Nout (1) Nband (2) parameter
338 else:
339 L = np.concatenate((-A2,-A3),1).transpose(2,0,1) # axis (0) Nout (1) Nband (2) parameter
341 return L,R
343 return LR,D,lsq_err,H,Npar
344def __parameterize_notrans(Nin,Nout,Nband,order):
345 """
346 Generate parameterization functions for estimation not including transient removal.
347 Parameterizations are passed into __solveLRM to generate local models
348 """
349 # Components of f, because f can be iterated
350 # Radius and exponent matrices. r is radius, s power, and A/B/D component matrices
351 N_half_band = (Nband-1)/2
352 r = np.arange(-N_half_band,N_half_band+1)
353 rsA = r**np.arange(order[0]+1)[:,np.newaxis] # axis. (0) power (1) radius
354 rsD = r**np.arange(order[2]+1)[:,np.newaxis]
356 # Define theta break points
357 bA = np.arange(Nout*Nin*(order[0]+1)) # transfer numerator
358 bD = np.arange(bA[-1]+1,bA[-1]+Nout*order[2]+1) # denominator (wrap Nout first, then s)
359 Npar = bA[-1]+Nout*order[2]+1
361 def lsq_err(U,Y,t):
362 t = t.reshape(-1,order='F') # transform from SO format to MO
364 # Build matrices
365 A = np.einsum('ij,jk->ik',t[bA].reshape(-1,order[0]+1,order='F'),rsA).reshape(Nout,Nin,-1,order='F')
366 B = np.zeros((Nout,Nband)) # axis (0) Nout (1) Nfreq
367 catmat = np.concatenate( (np.ones((Nout,1)),t[bD][:,np.newaxis]),0 )
368 D = np.einsum('ij,jk->ik',catmat.reshape(Nout,-1,order='F'),rsD) # axis. (1) Nout (2) Nfreq
370 # Compute sum of squared residuals
371 Ysim = (np.einsum('jkl,kl->jl',A,U) + B)/D
372 return (abs(Y-Ysim)**2).sum()
374 def H(t,r_inq):
375 t = t.reshape(-1,order='F') # transform from SO format to MO
377 # Radius and exponent matrices
378 rsA_inq = r_inq**np.arange(order[0]+1)[:,np.newaxis] # axis. (0) power (1) radius
379 rsD_inq = r_inq**np.arange(order[2]+1)[:,np.newaxis]
381 # Build matrices
382 A = np.einsum('ij,jk->ik',t[bA].reshape(-1,order[0]+1,order='F'),rsA_inq).reshape(Nout,Nin,-1,order='F') # axis. (1) Nout (2) Nin (3) Nfreq
383 catmat = np.concatenate( (np.ones((Nout,1)),t[bD][:,np.newaxis]),0 )
384 D = np.einsum('ij,jk->ik',catmat.reshape(Nout,-1,order='F'),rsD_inq) # axis. (1) Nout (2) Nfreq
386 return A/D[:,np.newaxis,:]
388 # Build linearized (SK) transfer matrix L(U,y)*t + R(U,Y) = Dp^(-1) (D(t)*Y - N(t)*U - M(t))
389 # Denominator matrix for SK iterations
390 def D(t):
391 return 1 + np.einsum( 'ij,lj->li',r[:,np.newaxis]**np.arange(1,order[2]+1) , t[:,np.arange(-order[2],0)] ) # axis (0) Nout (1) Nband
393 def LR(U,Y):
394 # Compute residual
395 R = -Y # axis (0) Nout (1) Nband
397 # Compute components
398 r1 = r[:,np.newaxis]
399 A1 = Y.T.reshape(-1,1,Nout,order='F')* (r1**np.arange(1,order[2]+1))[:,:,np.newaxis] # D(t)*Y. Axis. (1) r (2) s (3) Nout
400 A2 = (U.T[:,:,np.newaxis] * r1[:,:,np.newaxis] **np.arange(order[0]+1)).reshape(Nband,Nin*(order[0]+1),1,order='F') # A(t)*U. Axis. (1) r (2) Nin*s (3) Nout
402 # Assemble into array
403 if order[-1]>0:
404 catmat = np.ones((1,1,Y.shape[0])) # A2 and A3 don't depend on Nout
405 L = np.concatenate((-A2*catmat,A1),1).transpose(2,0,1) # axis (0) Nout (1) Nband (2) parameter
406 else:
407 L = -A2.transpose(2,0,1) # axis (0) Nout (1) Nband (2) parameter
409 return L,R
411 return LR,D,lsq_err,H,Npar
413#%% Frequency parser
414def __parse_bands(N_half_band,max_from_center,f_in,f_out):
415 """
416 Local modeling is not constrained to the input frequency vector. This function
417 processes input parameters and parses f_in into bands to model (bands).
418 In addition, it returns query points (query) at which to export each
419 local model
420 """
421 df_out = f_out[1]-f_out[0] # query point step size
422 df_in = f_in[1]-f_in[0] # input frequencies steps size
423 inq_per_band = mt.floor(2*max_from_center/df_out) + 1 # query points per band
425 # Assign export indices to band numbers
426 n_lo = (f_out <= f_in[N_half_band]+max_from_center).sum()
427 n_hi = (f_out >= f_in[-N_half_band-1]-max_from_center).sum()
428 n_mid = f_out.size - n_lo - n_hi
430 # Parsing is different if n_lo and n_hi > 0
431 is_lo = n_lo > 0
432 is_hi = n_hi > 0
434 # Get band assignment numbers for each entry in f_out
435 temp = np.ones((inq_per_band,1))*np.arange(is_lo,mt.ceil(n_mid/inq_per_band)+1)
436 temp = temp.reshape(-1,order='F')[:n_mid]
437 inq_assign = np.concatenate(( np.zeros(n_lo) , temp , (temp.max()+1)*np.ones(n_hi) ))
438 inq_assign = inq_assign.astype(dtype='int32')
440 # Get indices of center points and build bands
441 c_inds_hz = np.array([f_out[inq_assign==k].mean() for k in range(is_lo,inq_assign.max()+1-is_hi)]) # ideal center locations of each export
442 c_inds = ((c_inds_hz-f_in[0])/df_in).round() # actual center locations of each export on f_in
443 if abs(c_inds_hz - f_in[c_inds.astype('int32')]).max() > max_from_center:
444 warnings.warn("""f_in is too sparse for f_out specification, resulting in LM extrapolation
445 beyond the export_ratio. Consider using default specification or increasing bandwidth.""")
446 lo_ind = np.array([N_half_band]) if is_lo else np.ones(0)
447 hi_ind = np.array([f_in.size-1-N_half_band]) if is_hi else np.ones(0)
448 c_inds = np.concatenate(( lo_ind,c_inds,hi_ind )).astype(dtype='int32')
449 bands = np.arange(-N_half_band,N_half_band+1) + c_inds[:,np.newaxis] # axis (0) center ind (1) r
450 bands = bands.astype(dtype='int32')
452 # Get normalized extraction radii
453 query = [(f_out[inq_assign==k] - ind)/df_in for k,ind in enumerate(f_in[c_inds])]
455 return bands,query
457#%% General functions
458def __round_odd(x):
459 y = int(x)
460 return y + int(y%2==0)
461def __lstsq_scipy(A,b): # extension of np.linalg.lstsq to 3-D arrays. Not vectorized.
462 X = np.empty( np.array(A.shape)[[0,-1]] , dtype=complex ) # axis. (0) Nout [1] Npar
463 for k,Ak in enumerate(A):
464 X[k] = sp.linalg.lstsq(Ak,b[k])[0]
465 return X
466def __lstsq_numpy(A,b): # extension of np.linalg.lstsq to 3-D arrays. Not vectorized.
467 X = np.empty( np.array(A.shape)[[0,-1]] , dtype=complex ) # axis. (0) Nout [1] Np
468 for k,Ak in enumerate(A):
469 X[k] = np.linalg.lstsq(Ak,b[k],rcond=None)[0]
470 return X
472"""
473APPENDIX
474Note 1: Windows OS restricts parallel processes to 61
476Note 2: linalg.lstsq in scipy and numpy vary wildly in relative speed. In some cases,
477numpy is 10x faster. For other problem sizes, scipy is 20x faster. Both are tested,
478and the faster algorithm is selected.
479"""