deploy
In[]:=
Tue 8 Oct 2019 14:14:46

Cheap version of B3

SeedRandom[1];​​d=5;​​normalize[x_]:=x/Total[x,10];​​u=normalize[RandomReal[{0,1},d]];​​Unprotect[D];​​D=DiagonalMatrix[u];​​A=D-Outer[Times,u,u];​​cheap=DiagonalMatrix[Sqrt[u]]-Outer[Times,u,
u
];​​Assert[cheap.cheapA]
In[]:=

Compare against symmetric decomposition

SeedRandom[1];​​Unprotect[D];​​d=5;​​normalize[x_]:=x/Total[x,10];​​u=normalize[RandomReal[{0,1},d]];​​D=DiagonalMatrix[u];​​A=D+uu;​​​​​​uu=Outer[Times,u,u];(*uu'*)​​​​iD=DiagonalMatrix[1/u];(*
-1
D
*)​​x=
u.iD.u+1
-1
u.iD.u
;​​B2=
1/2
D
+x
1/4
iD
.uu.
1/4
iD
;​​B3=
1/2
D
+xuu.
1/2
iD
;​​​​​​Print["error exact=",Norm[MatrixPower[symsqrt[A],2]-A]]​​Print["error rank1=",Norm[B2-symsqrt[A]]]​​Print["error rank1 fixed=",Norm[B3-symsqrt[A]]]​​Print["error diag=",Norm[MatrixPower[D,1/2]-symsqrt[A]]]​​
In[]:=
error exact=3.99455×
-16
10
error rank1=0.00876574
error rank1 fixed=0.0414601
error diag=0.219797

Errors of two approximations

Clear[doit];​​doit[idx_]:=Module{},​​d=5;​​u=normalize[RandomReal[{0,1},d]];​​D=DiagonalMatrix[u];​​A=D+uu;​​​​​​uu=Outer[Times,u,u];(*uu'*)​​​​iD=DiagonalMatrix[1/u];(*
-1
D
*)​​x=
u.iD.u+1
-1
u.iD.u
;​​B2=
1/2
D
+x
1/4
iD
.uu.
1/4
iD
;​​B3=
1/2
D
+xuu.
1/2
iD
;​​truth=symsqrt[A];​​{Norm[B2-truth],​​Norm[B3-truth]}​​;​​vals=doit/@Range[1000000];​​Count[Thread[vals[[All,1]]<vals[[All,2]]],True]
In[]:=
897509
Out[]=
​