Switch to QR-based matrix inversion, bug fixes, bump version number

Project.toml

@@ -1,7 +1,7 @@
 name = "WildBootTests"
 uuid = "65c2e505-86ba-4c19-93f1-95506c1443d5"
 authors = ["droodman <[email protected]>"]
-version = "0.9.1"
+version = "0.9.2"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/WRE.jl

@@ -581,7 +581,8 @@ function HessianFixedkappa!(o::StrBootTest{T}, dest::AbstractMatrix{T}, is::Vect
 
 		if _jk
 			!iszero(κ) &&
-				(dest[row,1] = dot(i>0 ? o.Repl.XZ[:,i] : o.Repl.Xy₁par, j>0 ? o.Repl.V[:,j] : o.Repl.invXXXy₁par))
+				(dest[row,1] = dot(i>0 ? o.Repl.XZ[:,i] : o.Repl.Xy₁par, 
+				                   j>0 ? o.Repl.V[ :,j] : o.Repl.invXXXy₁par))
 			!isone(κ) &&
 				(dest[row,1] = iszero(κ) ? o.Repl.YY[i+1,j+1] : κ * dest[row,1] + (1 - κ) * o.Repl.YY[i+1,j+1])
 		end
@@ -770,7 +771,7 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T
 		if o.null
 			o.numerWRE .= view(o.β̈sAs,1:1,:) .+ (o.Repl.Rt₁ - o.r) / o.Repl.RRpar
 		else
-			o.numerWRE .=  view(o.β̈sAs,1:1,:) .- view(o.DGP.β̈ ,:,1)
+			o.numerWRE .= view(o.β̈sAs,1:1,:) .- view(o.DGP.β̈ ,:,1)
 			isone(w) && (o.numerWRE[1] = o.β̈sAs[1] + (o.Repl.Rt₁[1] - o.r[1]) / o.Repl.RRpar[1])
 		end
 
@@ -817,15 +818,15 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T
 			!iszero(o.fuller) && (o.κWRE .-= o.fuller / (o._Nobs - o.kX))
 
 			o.As .= o.κWRE .* view(o.YPXY✻, 2:o.Repl.kZ+1, :, 2:o.Repl.kZ+1) .+ (1 .- o.κWRE) .* view(o.YY✻, 2:o.Repl.kZ+1, :, 2:o.Repl.kZ+1)
-			invNaN!(o.As)
+			invsym!(o.As)
 			t✻!(view(o.β̈s,:,:,1:1), o.As, o.κWRE .* view(o.YPXY✻, 2:o.Repl.kZ+1, :, 1) .+ (1 .- o.κWRE) .* view(o.YY✻, 2:o.Repl.kZ+1, :,  1))
 		else
 			HessianFixedkappa!(o, o.δnumer, collect(1:o.Repl.kZ), 0, o.κ, _jk)
 			@inbounds for i ∈ 1:o.Repl.kZ
 				HessianFixedkappa!(o, view(o.As, 1:i, :, i), collect(1:i), i, o.κ, _jk)
 			end
 			symmetrize!(o.As)
-			invNaN!(o.As)
+			invsym!(o.As)
 			t✻!(view(o.β̈s,:,:,1:1), o.As, view(o.δnumer,:,:,1:1))
 		end
 
@@ -871,7 +872,7 @@ function MakeWREStats!(o::StrBootTest{T}, w::Integer) where T
 			if o.sqrt
 				@storeWtGrpResults!(o.dist, o.numerWRE ./ sqrtNaN.(dropdims(o.denomWRE; dims=3)))
 			else
-				invNaN!(o.denomWRE)
+				invsym!(o.denomWRE)
 				_numer = view(o.numerWRE,:,:,1:1)
 				@storeWtGrpResults!(o.dist, dropdims(_numer'o.denomWRE*_numer; dims=3))  # hand-code for 2-dimensional?  XXX allocations
 			end

src/WildBootTests.jl

@@ -95,6 +95,6 @@ function UpdateBootstrapcDenom!(o::StrBootTest{T} where T)
 	nothing
 end
 
-# include("precompiler.jl")
+include("precompiler.jl")
 
 end

src/estimators.jl

@@ -64,7 +64,7 @@ function InitVarsOLS!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstrac
 	end
 
 	o.invH = (pinv(H))
-  R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.invH : (o.R₁perp * invNaN(o.R₁perp'H*o.R₁perp) * o.R₁perp')  # for DGP regression
+  R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.invH : (o.R₁perp * invsym(o.R₁perp'H*o.R₁perp) * o.R₁perp')  # for DGP regression
 	o.β̈₀ = R₁AR₁ * X₁y₁
 	o.∂β̈∂r = R₁AR₁ * H * o.R₁invR₁R₁ - o.R₁invR₁R₁
 
@@ -79,17 +79,17 @@ function InitVarsOLS!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstrac
 				v = view(parent.X₁, S,:)
 				o.invMjk[g] = - v * R₁AR₁ * v'
 				o.invMjk[g][1:length(S)+1:length(S)^2] .+= one(T)  # add I
-				o.invMjk[g] .= invNaN(o.invMjk[g])
+				o.invMjk[g] .= invsym(o.invMjk[g])
 			end
 		else
 			!isdefined(o, :S✻XX) && (o.S✻XX = panelcross(parent.X₁, parent.X₁, parent.info✻))
 			_H = reshape(H, (parent.kX, 1, parent.kX)) .- o.S✻XX
-			_invH = iszero(nrows(o.R₁perp)) ? invNaN(_H) : o.R₁perp * invNaN(o.R₁perp' * _H * o.R₁perp) * o.R₁perp'
+			_invH = iszero(nrows(o.R₁perp)) ? invsym(_H) : o.R₁perp * invsym(o.R₁perp' * _H * o.R₁perp) * o.R₁perp'
 			o.XinvHjk = [view(parent.X₁, parent.info✻[g],:) * view(_invH,:,g,:) for g ∈ 1:parent.N✻]
 		end
 	end
 
-  o.A = iszero(nrows(Rperp)) ? o.invH : (Rperp * invNaN(Rperp'H*Rperp) * Rperp')  # for replication regression
+  o.A = iszero(nrows(Rperp)) ? o.invH : (Rperp * invsym(Rperp'H*Rperp) * Rperp')  # for replication regression
   o.AR = o.A * parent.R'
   (parent.scorebs || parent.robust) && (o.XAR = parent.X₁ * o.AR)
 	nothing
@@ -123,8 +123,8 @@ function InitVarsARubin!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
 	end
 
   H = ([X₁X₁ X₂X₁' ; X₂X₁ X₂X₂])
-  o.A = invNaN(H)
-  R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.A : (o.R₁perp * invNaN(o.R₁perp'H*o.R₁perp) * o.R₁perp')
+  o.A = invsym(H)
+  R₁AR₁ = iszero(nrows(o.R₁perp)) ? o.A : (o.R₁perp * invsym(o.R₁perp'H*o.R₁perp) * o.R₁perp')
 	o.β̈₀   = R₁AR₁ * [X₁y₁ ; X₂y₁]
 	o.∂β̈∂r = R₁AR₁ * [X₁Y₂ ; X₂Y₂]
 
@@ -145,12 +145,12 @@ function InitVarsARubin!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
 											v₂ * (@view negR₁AR₁[parent.kX₁+1:end, 1:parent.kX₁    ]) * v₁' +
 											v₂ * (@view negR₁AR₁[parent.kX₁+1:end, parent.kX₁+1:end]) * v₂'
 				o.invMjk[g][1:length(S)+1:length(S)^2] .+= one(T)  # add I
-				o.invMjk[g] .= invNaN(o.invMjk[g])
+				o.invMjk[g] .= invsym(o.invMjk[g])
 			end
 		else
 			!isdefined(o, :S✻XX) && (o.S✻XX = [[S✻X₁X₁ S✻X₂X₁'] ; [S✻X₂X₁ S✻X₂X₂]])
 			_H = reshape(H, (parent.kX, 1, parent.kX)) .- o.S✻XX
-			_invH = iszero(nrows(o.R₁perp)) ? invNaN(_H) : o.R₁perp * invNaN(o.R₁perp' * _H * o.R₁perp) * o.R₁perp'
+			_invH = iszero(nrows(o.R₁perp)) ? invsym(_H) : o.R₁perp * invsym(o.R₁perp' * _H * o.R₁perp) * o.R₁perp'
 			o.XinvHjk = [(S = parent.info✻[g]; X₁₂B(view(parent.X₁, S,:), view(parent.X₂, S,:), view(_invH,:,g,:))) for g ∈ 1:parent.N✻]
 		end
 	end
@@ -206,7 +206,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 
 		ZperpZperp, o.invZperpZperp, _invZperpZperp = invsymcrossjk(o.Zperp, parent.info✻)
 
-    _invZperpZperpZperpX₁  = _invZperpZperp * _ZperpX₁
+		_invZperpZperpZperpX₁  = _invZperpZperp * _ZperpX₁
     _invZperpZperpZperpX₂  = _invZperpZperp * _ZperpX₂
     _invZperpZperpZperpZ   = _invZperpZperp * _ZperpZ
     _invZperpZperpZperpZR₁ = _invZperpZperp * _ZperpZR₁
@@ -244,7 +244,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 
     o.XY₂ⱼₖ = [_X₁Y₂ ; _X₂Y₂]
     o.XXⱼₖ  = [_X₁X₁ ; _X₂X₁ ;;; _X₂X₁' ;  _X₂X₂]
-		o.invXXⱼₖ = invNaN(o.XXⱼₖ)
+		o.invXXⱼₖ = invsym(o.XXⱼₖ)
     o.H_2SLSⱼₖ = o.XZⱼₖ'o.invXXⱼₖ * o.XZⱼₖ
     (!isone(o.κ) || o.liml) && (o.H_2SLSmZZⱼₖ = o.H_2SLSⱼₖ - o.ZZⱼₖ)
 
@@ -276,7 +276,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 		end
 		ZperpZ   = o.Zperp'o.Zpar
 		o.restricted && (ZperpZR₁ = o.Zperp'o.ZR₁)
-		o.invZperpZperp = invNaN(o.Zperp'o.Zperp)
+		o.invZperpZperp = invsym(o.Zperp'o.Zperp)
 	end
 
 	if !o.isDGP && parent.WREnonARubin
@@ -307,7 +307,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 
 			X₂X₁ = o.X₂'o.X₁
 			o.XX = [o.X₁'o.X₁ X₂X₁' ; X₂X₁ o.X₂'o.X₂]
-			o.invXX = invNaN(o.XX)
+			o.invXX = invsym(o.XX)
 			X₁Y₂ = o.X₁'o.Y₂
 			X₂Y₂ = o.X₂'o.Y₂
 			o.XY₂ = [X₁Y₂ ; X₂Y₂]
@@ -362,7 +362,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 			S✻⋂X₂Zperp = panelcross(parent.X₂, o.Zperp, parent.info✻⋂)
 
 			if !(parent.jk && parent.WREnonARubin)
-				o.invZperpZperp = iszero(o.kZperp) ? Matrix{T}(undef,0,0) : invNaN(sumpanelcross(o.S✻⋂ZperpZperp))
+				o.invZperpZperp = iszero(o.kZperp) ? Matrix{T}(undef,0,0) : invsym(sumpanelcross(o.S✻⋂ZperpZperp))
 				ZperpX₁ = sumpanelcross(S✻⋂X₁Zperp)'
 				ZperpX₂ = sumpanelcross(S✻⋂X₂Zperp)'
 			end
@@ -385,7 +385,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 			X₁X₁ .-= ZperpX₁'o.invZperpZperp * ZperpX₁
 			X₂X₂ .-= ZperpX₂'o.invZperpZperp * ZperpX₂
 			o.XX = ([X₁X₁ X₂X₁' ; X₂X₁ X₂X₂])
-			o.invXX = invNaN(o.XX)
+			o.invXX = invsym(o.XX)
 			o.kX = ncols(o.XX)
 			o.kX₁ = ncols(X₂X₁)
 
@@ -419,6 +419,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 		  o.Y₂y₁ .-= ZperpY₂'o.invZperpZperpZperpy₁
 			o.Y₂Y₂ .-= ZperpY₂'o.invZperpZperpZperpY₂
 			o.X₂y₁ .-= ZperpX₂'o.invZperpZperpZperpy₁
+
 			o.X₁y₁ .-= ZperpX₁'o.invZperpZperpZperpy₁
 			o.y₁y₁  -= Zperpy₁'o.invZperpZperpZperpy₁
 
@@ -498,7 +499,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 	o.Zy₁par     = copy(o.Zy₁)
 	o.y₁pary₁par = copy(o.y₁y₁)
 	o.Xy₁par     = [o.X₁y₁ ; o.X₂y₁]
-	((parent.jk && !o.isDGP) || (!parent.jk && (parent.granular || parent.scorebs))) &&
+	((parent.jk && !o.isDGP) || (parent.granular || parent.scorebs)) &&
 		(o.y₁par   = copy(o.y₁))
 
 	o.V =  o.invXX * o.XZ  # in 2SLS case, estimator is (V' XZ)^-1 * (V'Xy₁). Also used in k-class and liml robust VCV by Stata convention
@@ -509,7 +510,7 @@ function InitVarsIV!(o::StrEstimator{T}, parent::StrBootTest{T}, Rperp::Abstract
 			if o.liml
 				o.invHⱼₖ = Array{T,3}(undef, o.kZ, parent.N✻, o.kZ)
 			else
-				o.invHⱼₖ = invNaN(isone(o.κ) ? o.H_2SLSⱼₖ : o.ZZⱼₖ + o.κ * o.H_2SLSmZZⱼₖ)
+				o.invHⱼₖ = invsym(isone(o.κ) ? o.H_2SLSⱼₖ : o.ZZⱼₖ + o.κ * o.H_2SLSmZZⱼₖ)
 			end
 		end
 
@@ -554,10 +555,10 @@ end
 
 function MakeH!(o::StrEstimator{T}, parent::StrBootTest{T}, makeXAR::Bool=false) where T
   H = isone(o.κ) ? o.H_2SLS : o.ZZ + o.κ * o.H_2SLSmZZ
-  o.invH = invNaN(H)
+  o.invH = invsym(H)
 
   if makeXAR  # for replication regression in score bootstrap of IV/GMM
-	  o.A = ncols(o.Rperp)>0 ? (o.Rperp * invNaN(o.Rperp'H*o.Rperp) * o.Rperp') : o.invH
+	  o.A = ncols(o.Rperp)>0 ? (o.Rperp * invsym(o.Rperp'H*o.Rperp) * o.Rperp') : o.invH
 	  o.AR = o.A * (o.Rpar'parent.R')
 	  o.XAR = X₁₂B(o.X₁, o.X₂, o.V * o.AR)
   end
@@ -598,7 +599,7 @@ function EstimateIV!(o::StrEstimator{T}, parent::StrBootTest{T}, _jk::Bool, r₁
 
 		if o.liml
 			o.YPXY = ([[o.invXXXy₁par'o.Xy₁par] o.ZXinvXXXy₁par' ; o.ZXinvXXXy₁par  o.H_2SLS])
-	    o.κ = 1/(1 - real(eigvalsNaN(invNaN(o.YY) * o.YPXY)[1]))  # like Fast & Wild (81), but more stable, at least in Mata
+	    o.κ = 1/(1 - real(eigvalsNaN(invsym(o.YY) * o.YPXY)[1]))  # like Fast & Wild (81), but more stable, at least in Mata
 	    !iszero(o.fuller) && (o.κ -= o.fuller / (parent._Nobs - parent.kX))
 	    MakeH!(o, parent)
 		end
@@ -611,9 +612,9 @@ function EstimateIV!(o::StrEstimator{T}, parent::StrBootTest{T}, _jk::Bool, r₁
 
 			if o.liml
 				o.YPXYⱼₖ .= [o.invXXXy₁parⱼₖ'o.Xy₁parⱼₖ ; o.ZXinvXXXy₁parⱼₖ ;;; o.ZXinvXXXy₁parⱼₖ' ; o.H_2SLSⱼₖ]
-				o.κⱼₖ .= reshape(one(T) ./ (one(T) .- getindex.(real.(eigvalsNaN.(each(invNaN(o.YYⱼₖ) * o.YPXYⱼₖ))), 1)), (1,:,1))
+				o.κⱼₖ .= reshape(one(T) ./ (one(T) .- getindex.(real.(eigvalsNaN.(each(invsym(o.YYⱼₖ) * o.YPXYⱼₖ))), 1)), (1,:,1))
 				!iszero(o.fuller) && (o.κⱼₖ .-= reshape(o.fuller ./ (o.Nobsⱼₖ .- parent.kX)), (1,:,1))
-				o.invHⱼₖ .= o.ZZⱼₖ .+ o.κⱼₖ .* o.H_2SLSmZZⱼₖ; invNaN!(o.invHⱼₖ)
+				o.invHⱼₖ .= o.ZZⱼₖ .+ o.κⱼₖ .* o.H_2SLSmZZⱼₖ; invsym!(o.invHⱼₖ)
 				o.β̈ⱼₖ .= o.κⱼₖ .* (o.ZXinvXXXy₁parⱼₖ .- o.Zy₁parⱼₖ) .+ o.Zy₁parⱼₖ
 			else
 				if isone(o.κ)
@@ -680,7 +681,7 @@ function MakeResidualsIV!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
 
 	  Xu = o.Xy₁par - o.XZ * o.β̈  # after DGP regression, compute Y₂ residuals by regressing Y₂ on X while controlling for y₁ residuals, done through FWL
 	  negXuinvuu = Xu / -uu
-	  o.Π̈ = invNaN(o.XX + negXuinvuu * Xu') * (negXuinvuu * (o.Y₂y₁par - o.ZY₂'o.β̈ )' + o.XY₂)
+	  o.Π̈ = invsym(o.XX + negXuinvuu * Xu') * (negXuinvuu * (o.Y₂y₁par - o.ZY₂'o.β̈ )' + o.XY₂)
 		o.γ̈ = o.Rpar * o.β̈  + o.t₁ - parent.Repl.t₁
     o.Ü₂Ü₂ = o.Y₂Y₂ - (o.Π̈)'o.XY₂ - o.XY₂'o.Π̈ + (o.Π̈)'o.XX*o.Π̈
 
@@ -710,7 +711,7 @@ function MakeResidualsIV!(o::StrEstimator{T}, parent::StrBootTest{T}) where T
 
 				Xu .= view(o.Xy₁parⱼₖ,:,g,:); t✻minus!(Xu, view(o.XZⱼₖ,:,g,:), view(o.β̈ⱼₖ,:,g,1))  # after DGP regression, compute Y₂ residuals by regressing Y₂ on X while controlling for y₁ residuals, done through FWL
 				negXuinvuu .= Xu ./ -uu
-				Π̈ⱼₖ .= invNaN(view(o.XXⱼₖ,:,g,:) + negXuinvuu * Xu') * (negXuinvuu * (view(o.Y₂y₁parⱼₖ,:,g,:) - view(o.ZY₂ⱼₖ,:,g,:)'view(o.β̈ⱼₖ,:,g,1))' + view(o.XY₂ⱼₖ,:,g,:))
+				Π̈ⱼₖ .= invsym(view(o.XXⱼₖ,:,g,:) + negXuinvuu * Xu') * (negXuinvuu * (view(o.Y₂y₁parⱼₖ,:,g,:) - view(o.ZY₂ⱼₖ,:,g,:)'view(o.β̈ⱼₖ,:,g,1))' + view(o.XY₂ⱼₖ,:,g,:))
 				X₁₂Bminus!(view(o.Ü₂,S,:), view(o.X₁ⱼₖ,S,:), view(o.X₂ⱼₖ,S,:), Π̈ⱼₖ)
 				t✻plus!(view(o.u⃛₁,S), view(o.Ü₂,S,:), o.RparY * view(o.β̈ⱼₖ,:,g,1) + _t₁Y)
 			end

src/interface.jl

@@ -14,7 +14,7 @@ struct BootTestResult{T}
   b::Vector{T}
   V::Matrix{T}
   auxweights::Union{Nothing,Matrix{T}}
-  o::StrBootTest
+  # o::StrBootTest
 end
 
 "Return test statistic"
@@ -155,7 +155,7 @@ function __wildboottest(
 	                  o.p, padj, o.B, o.BFeas, o.N✻, o.dof, o.dof_r, plot, peak, ci,
 	                  getdist(o, diststat),
 	                  getb(o), getV(o),
-	                  getauxweights && reps>0 ? getv(o) : nothing , o)
+	                  getauxweights && reps>0 ? getv(o) : nothing #=, o=#)
 end
 
 vecconvert(T::DataType, X) = isa(X, Vector{T}) ? X : Vector{T}(reshape(X, size(X,1)))

src/nonWRE.jl

@@ -289,7 +289,7 @@ function MakeNonWREStats!(o::StrBootTest{T}, w::Integer) where T
 			@storeWtGrpResults!(o.dist, o.numerw ./ sqrtNaN.(o.denom[1,1]))
 			isone(w) &&
 				(o.statDenom = hcat(o.denom[1,1][1]))  # original-sample denominator
-		elseif o.dof==2  # hand-code 2D numer'invNaN(denom)*numer
+		elseif o.dof==2  # hand-code 2D numer'invsym(denom)*numer
 			t1 = view(o.numerw,1,:)'; t2 = view(o.numerw,2,:)'; t12 = t1.*t2
 			@storeWtGrpResults!(o.dist, (t1.^2 .* o.denom[2,2] .- 2 .* t12 .* o.denom[2,1] .+ t2.^2 .* o.denom[1,1]) ./ (o.denom[1,1].*o.denom[2,2] .- o.denom[2,1].^2))
 			isone(w) &&

src/structs.jl

@@ -54,22 +54,22 @@ end
 
 function *(X::AbstractMatrix{T}, Y::DesignerMatrix{T}) where T
 	if Y.type==regular
-		view(Base.:*(X, Y.M),:,collect(1:size(Y.M,2)))  # collect() makes the column slicer Vector{Int64}, for type consistency
+		Base.:*(X, Y.M)  # collect() makes the column slicer Vector{Int64}, for type consistency
 	elseif Y.type==identity
-		view(X,:,collect(1:size(X,2)))
+		X
 	elseif length(Y.p)==Y.size[2]
-		 view(X,:,Y.p[Y.q])
+		 X[:,Y.p[Y.q]]
 	else
 		dest = zeros(T, size(X,1), Y.size[2])  # Y is selection matrix with some cols all 0
 		if length(dest)>0
-			@inbounds for j ∈ axes(Y.p,1)
+			@inbounds for j ∈ eachindex(axes(Y.p,1))
 				pⱼ, qⱼ = Y.p[j], Y.q[j] 
 				@tturbo warn_check_args=false for i ∈ indices((dest,X),1)
-				 dest[i,qⱼ] = X[i,pⱼ]
+					dest[i,qⱼ] = X[i,pⱼ]
 				end
 			end
 		end
-		view(dest,:,collect(1:size(dest,2)))
+		dest
 	end
 end
 function *(X::AbstractArray{T,3}, Y::DesignerMatrix{T}) where T
@@ -170,7 +170,7 @@ mutable struct StrEstimator{T<:AbstractFloat}
 	ZparX::Matrix{T}
 	invZperpZperpZperpX₁::Matrix{T}; invZperpZperpZperpX₂::Matrix{T}; invZperpZperpZperpy₁::Vector{T}; invZperpZperpZperpY₂::Matrix{T}; S✻UY₂::Matrix{T}; invZperpZperpZperpZpar::Matrix{T}; invZperpZperpZperpZR₁::Matrix{T}
 	Ü₂Ü₂::Matrix{T}; γ̈X::Vector{T}; γ̈Y::Vector{T}; γ⃛::Vector{T}; Xȳ₁::Vector{T}; ȳ₁ȳ₁::T; XÜ₂::Matrix{T}; ȳ₁Ü₂::Matrix{T}; Ȳ₂::Matrix{T}; ȳ₁::Vector{T}
-	Xpar₁toZparX::DesignerMatrix{T}; X₁noFWL::DesignerProduct{T}
+	Xpar₁toZparX::DesignerMatrix{T}; X₁noFWL::Matrix{T} #=DesignerProduct{T}=#
 
 	X₁ⱼₖ::Matrix{T}; X₂ⱼₖ::Matrix{T}; y₁ⱼₖ::Vector{T}; Y₂ⱼₖ::Matrix{T}; Zⱼₖ::Matrix{T}; ZR₁ⱼₖ::Matrix{T}; Y₂y₁ⱼₖ::Array{T,3}; X₂y₁ⱼₖ::Array{T,3}; X₁y₁ⱼₖ::Array{T,3}; Zy₁ⱼₖ::Array{T,3}; XZⱼₖ::Array{T,3}; ZZⱼₖ::Array{T,3}; ZY₂ⱼₖ::Array{T,3}; y₁y₁ⱼₖ::Array{T,3}; XY₂ⱼₖ::Array{T,3}; invXXⱼₖ::Array{T,3}; XXⱼₖ::Array{T,3}; X₁ZR₁ⱼₖ::Array{T,3}; X₂ZR₁ⱼₖ::Array{T,3}; ZZR₁ⱼₖ::Array{T,3}; twoZR₁y₁ⱼₖ::Array{T,3}; ZR₁ZR₁ⱼₖ::Array{T,3}; ZR₁Y₂ⱼₖ::Array{T,3} 
 	Y₂y₁parⱼₖ::Array{T,3}; Zy₁parⱼₖ::Array{T,3}; y₁pary₁parⱼₖ::Array{T,3};	Xy₁parⱼₖ::Array{T,3}; y₁parⱼₖ::Vector{T}; H_2SLSⱼₖ::Array{T,3}; H_2SLSmZZⱼₖ::Array{T,3}; invHⱼₖ::Array{T,3}

src/utilities.jl

@@ -20,22 +20,7 @@ end
 # iszero(nrows(X)) && (return Symmetric(X))
 # X, ipiv, info = LinearAlgebra.LAPACK.sytrf!('U', Matrix(X))
 # iszero(info) && LinearAlgebra.LAPACK.sytri!('U', X, ipiv)
-function invNaN(X)
-	iszero(length(X)) && (return X)
-	try
-		pinv(X)
-	catch _
-		fill(eltype(X)(NaN), size(X))
-	end
-end
-function invsym(X)
-	iszero(length(X)) && (return X)
-	try
-		pinv(Symmetric(X))
-	catch _
-		fill(eltype(X)(NaN), size(X))
-	end
-end
+@inline invsym(X) = X \ I
 
 eigvalsNaN(X) =
 	try
@@ -845,16 +830,16 @@ function t✻minus!(dest::AbstractArray{T,3}, A::AbstractVecOrMat{T}, B::Abstrac
 end
 
 # in-place inverse of a set of symmetric matrices
-function invNaN!(A::Array{T,3}) where T
+function invsym!(A::Array{T,3}) where T
 	@inbounds for g ∈ eachindex(axes(A,2))
-		A[:,g,:] = invNaN(@view A[:,g,:])
+		A[:,g,:] = invsym(@view A[:,g,:])
 	end
 	nothing
 end
-function invNaN(A::Array{T,3}) where T
+function invsym(A::Array{T,3}) where T
 	dest = similar(A)
 	@inbounds for g ∈ eachindex(axes(A,2))
-		dest[:,g,:] = invNaN(@view A[:,g,:])
+		dest[:,g,:] = invsym(@view A[:,g,:])
 	end
 	dest
 end
@@ -873,19 +858,19 @@ function crossjk(A::VecOrMat{T}, B::Vector{T}, info::Vector{UnitRange{Int64}}) w
 	(vec(sumt), t)
 end
 
-# given data matrix X and cluster-defining info vector, compute X'X, invNaN(X'X), and delete-g invNaN(X'X)'s efficiently
+# given data matrix X and cluster-defining info vector, compute X'X, invsym(X'X), and delete-g invsym(X'X)'s efficiently
 function invsymcrossjk(X::Matrix{T}, info::Vector{UnitRange{Int64}}) where T
 	SXX =  panelcross(X,X,info)
 	XX = sumpanelcross(SXX)
-	invXX = invNaN(XX)
+	invXX = invsym(XX)
 	for (g,S) ∈ enumerate(info)
 		Xg = view(X,S,:)
 		if size(Xg,1) > size(Xg,2)
-			SXX[:,g,:] = invNaN(XX - Xg'Xg)
+			SXX[:,g,:] = invsym(XX - Xg'Xg)
 		else
 			XginvXX = Xg * invXX
 			tmp = XginvXX * Xg'; tmp -= I
-			SXX[:,g,:] = invXX; t✻minus!(view(SXX,:,g,:), XginvXX'invNaN(tmp), XginvXX)
+			SXX[:,g,:] = invXX; t✻minus!(view(SXX,:,g,:), XginvXX'invsym(tmp), XginvXX)
 		end
 	end
 	(XX, invXX, SXX)