In[]:=
Mon 4 Dec 2023 18:35:55
Util
Util
In[]:=
Out[]=
AssertNorm[sylvesterSpectral[a,b,c]-LyapunovSolve[a,b,c]]<
1
10
10
Expectations framework
Expectations framework
In[]:=
On[Assert];Unprotect[Ex];Clear[x,y,X,xx,yy,xx1,xx2];(*symmetricMatrix[d_,sym_]:=(Clear[sym];Array[sym[Min[#1,#2],Max[#1,#2]]&,{d,d}]);*)(*ChangedNov10*)symmetricMatrix[d_,sym_]:=(Clear[sym];Array[sym[Max[#1,#2],Min[#1,#2]]&,{d,d}]);arbitaryMatrix[d_,sym_]:=(Clear[sym];Array[sym,{d,d}]);subVars[p_]:=With[{vars=Variables[p]},Thread[varsRandomInteger[{-5,5},Length[vars]]]]Ex[expr_]:=Expectation[expr,xdist];(*Expectationwithrespecttox*)Exy[expr_]:=Expectation[expr,{xdist,ydist}];(*Expectationwithrespecttox*)EX[expr_]:=Expectation[expr,{x1dist,x2dist}];Protect[Ex];(*Returnsx,y,A,As,H,Hs,a,as*)computeDerived:=Module[{},Clear[aa,hh];Clear[xx];Assert[Head[dist]~endsWith~"Distribution"];d=Length@RandomVariate[dist];n=d;(*symbolicvariablesforvectorandmatrixX*);x=Array[xx,d];y=Array[yy,d];x1=Array[xx1,d];x2=Array[xx2,d];(*batchofXs*);X={x1,x2};A=arbitaryMatrix[d,aa];As=symmetricMatrix[d,aa];symmetrizer=getSymmetrizerMatrix[n];H=arbitaryMatrix[d,hh];Hs=symmetricMatrix[d,hh];a=vec[A];as=vec[As];asVars=Flatten[As];ii=IdentityMatrix[d];ones=ConstantArray[1,d];vars=Variables[A]~Join~Variables[H];mu=Ex[x];X2=Ex[x⊗x];H=X2;X2X2=Ex[(x⊗x).(x⊗x)];X4=Ex[Outer[Times,x,x,x,x]];X4flat=Flatten[X4,{{1,2},{3,4}}];Assert[X4flatEx[(x⊗x)⊗(x⊗x)]];];(*Setupsymbolicandnumericmatricesforchecks*)subst:=(SeedRandom[1];Thread[varsRandomReal[{-1,1},{Length[vars]}]]);substInt:=(SeedRandom[1];Thread[varsRandomInteger[{-4,4},{Length[vars]}]]);(*setupdistribution*)setupDistEmpirical:=(d=2;n=d;X=RandomInteger[{-5,5},{d,20}];dist=EmpiricalDistribution[Transpose@X];computeDerived;);setupDistDataset[dataset_]:=Module[{},d=Last@Dimensions@dataset;n=d;X=dataset;dist=EmpiricalDistribution[X];computeDerived;];setupDistGaussian:=(d=3;n=d;Σ=RandomInteger[{-5,5},{d,d}];Σ=Σ.Transpose[Σ]+IdentityMatrix[d];μ=RandomInteger[{-5,5},{d}];Assert[Min[N@Eigenvalues[Σ]]>0];dist=MultinormalDistribution[μ,Σ];computeDerived);setupDistCenteredGaussian:=(d=3;n=d;Σ=RandomInteger[{-5,5},{d,d}];Σ=Σ.Transpose[Σ]+IdentityMatrix[d];μ=ConstantArray[0,{d}];Assert[Min[N@Eigenvalues[Σ]]>0];dist=MultinormalDistribution[Σ];computeDerived);setupDistDiagonalGaussian:=(d=4;n=d;Σ=DiagonalMatrix[Table[1/i,{i,1,d}]];μ=Table[0,d];Assert[Min[N@Eigenvalues[Σ]]>0];dist=MultinormalDistribution[Σ];computeDerived)
Added from linear-estimation-headers
Added from linear-estimation-headers
In[]:=
(******* Python interface *******)RegisterExternalEvaluator["Python","/Users/yaroslavvb/miniconda3/envs/kac/bin/python"];session=StartExternalSession["Python"];ExternalEvaluate[session,"import numpy as np"];saveDataCode="def f(data, fn): dataRoot='/Users/yaroslavvb/Library/CloudStorage/Dropbox/git0/kaczmarz/data' np.save(dataRoot+'/'+fn, data)";loadDataCode="def f(fn): dataRoot='/Users/yaroslavvb/Library/CloudStorage/Dropbox/git0/kaczmarz/data' return np.load(dataRoot+'/'+fn)";saveFunc=ExternalFunction[session,saveDataCode];loadFunc=ExternalFunction[session,loadDataCode];(* flips vector to make first coordinate positive *)signNormalize[vec_]:=vec*Sign[Total[vec]];topEvec[mat_]:=signNormalize@First@Eigenvectors[mat,1];topEvec[mat1_,mat2_]:=signNormalize@First@Eigenvectors[{mat1,mat2},1];topEval[mat1_]:=First@Eigenvalues[mat1,1];bottomEval[mat1_]:=Last@Eigenvalues[mat1];topEval[mat1_,mat2_]:=First@Eigenvalues[{mat1,mat2},1];bottomEval[mat1_,mat2_]:=Last@Eigenvalues[{mat1,mat2}];norm2[stuff_]:=With[{x=Flatten[stuff]},Total[x*x]]getCov[X_]:=With[{m=Length[X]},X.X/m];getCov4[X_]:=With[{norms=norm2/@X},(X).(norms*X)/Length[norms]];(* Symmetric square root *)symsqrt[m_]:=Module[{U,S,W},{U,S,W}=SingularValueDecomposition[m];U.Sqrt[S].Transpose[W]];(* Symmetric inverse square root *)isymsqrt[mat_]:=Module{u,e},{e,u}=Eigensystem[mat];u.DiagonalMatrix1,#]&,s];U.DiagonalMatrix[sdiv].Transpose[W]];norm2[stuff_]:=With[{x=Flatten[stuff]},Total[x x]];NN[expr_]:=N[expr,10]
e
.u;(* Symmetric inverse square root, using SVD for stability *)isymsqrtStable[m_]:=Module[{U,S,W,s,cutoff,sdiv},{U,S,W}=SingularValueDecomposition[m];s=Diagonal[S];(* float32 precision cutoff, borrowed from scipy pinv2 *);cutoff=Max[s]*Length[s]*$MachineEpsilon;sdiv=Map[If[#>cutoff,-1/2
#
General distribution formulas
General distribution formulas
Gaussian batch fourth moment formulas
Gaussian batch fourth moment formulas
Checking formula produced in https://stats.stackexchange.com/questions/589669/gaussian-fourth-moment-formulas
(unfishished) Follow-up is in https://stats.stackexchange.com/questions/589849/value-of-ext-a-x-b-xt-cx-for-gaussian-matrices-x
(unfinished) https://mathematica.stackexchange.com/questions/273893/automatically-deriving-formulas-for-wishart-moments
(unfishished) Follow-up is in https://stats.stackexchange.com/questions/589849/value-of-ext-a-x-b-xt-cx-for-gaussian-matrices-x
(unfinished) https://mathematica.stackexchange.com/questions/273893/automatically-deriving-formulas-for-wishart-moments
Dot product formulas
Dot product formulas
Normalized SNR formulas
Normalized SNR formulas
Gaussian Formulas
Gaussian Formulas
Matrix Formulas
Matrix Formulas
Gaussian Power iteration
Gaussian Power iteration
Gaussian Radius
Gaussian Radius
Upper bound from optimization
Upper bound from optimization
Matrix form
Matrix form
Improved bound
Improved bound
Isserlis factorization
Isserlis factorization
Cleaning notebook
Largest Files
Largest Files
Latest modified cells
Latest modified cells