Jan 29th.  Not sure what is below!

Needs["Quaternions`"]
HJtensorProduct[x_,y_]:= Module[
{QZero=Quaternion[0,0,0,0],Q,QTmp,X=Apply[Quaternion,x],Y=Apply[Quaternion,y]},
Table[
Q=QZero; Q[[i]]=1.0;
QTmp=QZero; QTmp[[i]]=1.0;
QTmp = X**QTmp;
QTmp=QTmp**Conjugate[Y];
Apply[List,QTmp],
{i,1,4}
]
]
HamiltonJacobi[S_,SVDFlag_]:=Module[
{T=ConstantArray[0,{3,3}]},
T[[1,1]]=0.25*(S[[1,1]]+S[[2,2]]-S[[3,3]]-S[[4,4]]);
T[[2,1]]=0.5*(S[[2,3]]+S[[1,4]]);
T[[3,1]]=0.5*(S[[2,4]]-S[[1,3]]);
T[[1,2]]=0.5*(S[[2,3]]-S[[1,4]]);
T[[2,2]]=0.25*(S[[1,1]]-S[[2,2]]+S[[3,3]]-S[[4,4]]);
T[[3,2]]=0.5*(S[[3,4]]+S[[1,2]]);
T[[1,3]]=0.5*(S[[2,4]]+S[[1,3]]);
T[[2,3]]=0.5*(S[[3,4]]-S[[1,2]]);
T[[3,3]]=0.25*(S[[1,1]]-S[[2,2]]-S[[3,3]]+S[[4,4]]);
TOut=T;
Switch[SVDFlag,
∞,{u,σ,v}=SingularValueDecomposition[T,1];{u,σ,v}=Map[Flatten,{u,σ,v}],
_Integer,{u,σ,v}= CheapDominantSVD[T, SVDFlag],
_, Print["InnappropriateSVDFlag ", SVDFlag]
];
u=Prepend[u,0];
v=Prepend[v,0];
u[[2]]+=1;
v[[2]]+=1;
u=u/Norm[u];
v=v/Norm[v];
R=HJtensorProduct[u,v]
]
CheapDominantSVD[A_,n_]:= Module[
{u,v,σ=0,B=Transpose[A].A},
v=A[[1]];
Do[
v=B.v;
σ=Norm[v];
v=v/σ,
{n}];
u=A.v;
u=u/Norm[u];
{u,(Sqrt[σ]),v}
]
CheapDominantShiftSVD[A_,n_]:= Module[
{u,v,σ=0,B=Transpose[A].A},
v=A[[All,1]];
Do[
v=(B+σ IdentityMatrix[3]).v;
σ=Norm[v];
v=v/σ,
{n}];
u=A.v;
u=u/Norm[u];
{u,(Sqrt[σ]),v}
]

Testing on single matrices with the exact eigenvalues.

A=RandomReal[{-1,1},{4,4}];
A=A+Transpose[A];
Q= HamiltonJacobi[A,∞];
Norm[IdentityMatrix[4]- Q.Transpose[Q]]
MatrixForm[Chop[Q.A.Transpose[Q]]]

Testing on bits of a bigger matrix.

n=8;
A=RandomReal[{-1,1},{n,n}];
A=A+Transpose[A];
ANew=A;
Indices =Sort[RandomSample[Range[1,n],4]];
Q= HamiltonJacobi[A[[Indices,Indices]],∞];
BigQ=SparseArray[Band[{1,1}]->1.0,{n,n}];
Do[
BigQ[[Indices[[i]],Indices[[j]]]]=Q[[i,j]],
{i, 4},{j,4}];
Norm[IdentityMatrix[n]- BigQ.Transpose[BigQ]]
ANew=BigQ.A.Transpose[BigQ];
MatrixForm[Chop[ANew]]

Trying to get organized.  If the index sets are disjoint the computations can be done in parallel

n=8;
A=RandomReal[{-1,1},{n,n}];
AOriginal=A+Transpose[A];
A=AOriginal;
Indices ={1,2,3,4};
Q= HamiltonJacobi[A[[Indices,Indices]],∞];
BigQ=SparseArray[Band[{1,1}]->1.0,{n,n}];
Do[
BigQ[[Indices[[i]],Indices[[j]]]]=Q[[i,j]],
{i, 4},{j,4}];
A=BigQ.A.Transpose[BigQ];
MatrixForm[Chop[A]]

Indices ={5,6,7,8};
Q= HamiltonJacobi[A[[Indices,Indices]],∞];
BigQ=SparseArray[Band[{1,1}]->1.0,{n,n}];
Do[
BigQ[[Indices[[i]],Indices[[j]]]]=Q[[i,j]],
{i, 4},{j,4}];
A=BigQ.A.Transpose[BigQ];
MatrixForm[Chop[A]]

It is not clear what to do next.  I guess I should see what happens when we do Hessenberg

n=12;
A=RandomReal[{-1,1},{n,n}];
AOriginal=A+Transpose[A];
A=AOriginal;
HJFlush[A_,Indices_]:= Module[
{Q,BigQ,n=Length[A]},
Q= HamiltonJacobi[A[[Indices,Indices]],∞];
BigQ=SparseArray[Band[{1,1}]->1.0,{n,n}];
Do[
BigQ[[Indices[[i]],Indices[[j]]]]=Q[[i,j]],
{i, 4},{j,4}];
Chop[BigQ.A.Transpose[BigQ]]
]
A=HJFlush[A,{1,2,3,4}];
A=HJFlush[A,{5,6,7,8}];
A=HJFlush[A,{5,6,7,8}];
A=HJFlush[A,{9,10,11,12}];
MatrixPlot[A]
Norm[Eigenvalues[A]-Eigenvalues[AOriginal]]

Try sweeping on a four diagonal matrix.

n=12;
A=SparseArray[{
Band[{1,1}]:> RandomReal[{-1,1}],
Band[{1,2}]:> RandomReal[{-1,1}],
Band[{1,3}]:>RandomReal[{-1,1}],
Band[{1,4}]:>RandomReal[{-1,1}],
Band[{2,1}]:> RandomReal[{-1,1}],
Band[{3,1}]:>RandomReal[{-1,1}],
Band[{4,1}]:>RandomReal[{-1,1}]
},
{n,n}];
AOriginal=A+Transpose[A];
A=AOriginal;
HJFlush[A_,Indices_]:= Module[
{Q,BigQ,n=Length[A]},
Q= HamiltonJacobi[A[[Indices,Indices]],∞];
BigQ=SparseArray[Band[{1,1}]->1.0,{n,n}];
Do[
BigQ[[Indices[[i]],Indices[[j]]]]=Q[[i,j]],
{i, 4},{j,4}];
Chop[BigQ.A.Transpose[BigQ]]
]
A=HJFlush[A,{1,2,3,4}];
A=HJFlush[A,{2,3,4,5}];
MatrixPlot[A]

ListPlot[{
Sort[Normal[Diagonal[AOriginal]]],
Sort[Normal[Diagonal[A]]],
Sort[Eigenvalues[A]]},
PlotStyle->{Red,Green,Blue}]

ListPlot3D[Abs[A]]

Transcribed From Jizhou's Matlab

Jizhou's SVD is quite a lot of arithmetic. Maybe I can beat it!

Fast3x3SVD1[A_, iter_]:= Module[
{i=0,B=Transpose[A].A,u,v=A[[All,1]]},
Do[
v=B.v;
v=(v/(Norm[v])),
{iter}];
u=A.v;
u=u/Norm[u];
{u,v}
]

Data=Table[
A=RandomReal[{-1,1},{3,3}];
{U,W,V}=SingularValueDecomposition[A];
{uM,vM}={U[[All,1]],V[[All,1]]};
{u2,v2}=Fast3x3SVD1[A,2];
{u4,v4}=Fast3x3SVD1[A,4];
{
Chop[ArcCos[((u2.uM)/(Norm[u2]Norm[uM]))]/π,(10^(-12))],
Chop[ArcCos[((v2.vM)/(Norm[v2]Norm[vM]))]/π,(10^(-12))],
Chop[ArcCos[((u4.uM)/(Norm[u4]Norm[uM]))]/π,(10^(-12))],
Chop[ArcCos[((v4.vM)/(Norm[v4]Norm[vM]))]/π,(10^(-12))]
},
{1000}];

ListPlot[
Table[
Sort[Data[[All,i]]],
{i,1,4}],
PlotStyle->{Red, Blue,Green,Pink}
]

A=RandomReal[{-1,1},{3,3}];
{U,W,V}=SingularValueDecomposition[A];
Sqrt[Diagonal[W]]
{eVal,eVec}=Eigensystem[A]

1+1

{U,W,V}=SingularValueDecomposition[A];
{eVal,eVec}=Eigensystem[A];
Sqrt[Diagonal[W]]
eVal

xx

Di

0.1*180/π

W

MatrixForm[U]

?Qua

Needs["Quaternions`"]
x=Quaternion[x1,x2,x3,x4];
y=Quaternion[y1,y2,y3,y4];
x+y
x**y
x/(4.0+ I)

?Quaternion

Break

TOut//MatrixForm

Solve[(x^3)+ a (x^2)+b x + c==0,x]//Simplify

TOut.Transpose[TOut]//MatrixForm

AA=Flatte

(B=A.Transpose[A])//MatrixForm

I am going to use
	(({{(Subscript[d, 1]), (Subscript[γ, 12]), (Subscript[γ, 13])}, {(Subscript[γ, 12]), (Subscript[d, 2]), (Subscript[γ, 23])}, {(Subscript[γ, 13]), (Subscript[γ, 23]), (Subscript[d, 3])}}))

I have a sym + semi def matrix. I want the largest eigenvector. I would love yto get this down to an involved
function of 3 scaled variables.  I have 
	(Subscript[d, i])≥0 
and
	(Subscript[d, i])+(Subscript[d, j])≥2| (Subscript[τ, ij])|
and 
	(Subscript[d, 1])+(Subscript[d, 2])+(Subscript[d, 3])≥2(| (Subscript[τ, 12])|+| (Subscript[τ, 13])|+| (Subscript[τ, 23])|)
The largest eigevalue and eigenvector are real. If I scale by 1/((Subscript[d, 1])+(Subscript[d, 2])+(Subscript[d, 3]))
I get a matrix with diagonal sum 1 and so
	0≤(Subscript[d, i])≤1 with (Subscript[d, 3])=1-(Subscript[d, 1])-(Subscript[d, 2])
	((|(Subscript[τ, ij])|)/((Subscript[d, 1])+(Subscript[d, 2])+(Subscript[d, 3])))≤(((Subscript[d, i])+(Subscript[d, j]))/(2((Subscript[d, 1])+(Subscript[d, 2])+(Subscript[d, 3]))))≤1/2
In other words the matrix has only 5 parameters and they are all of bounded range. My only remaining problem is that I need to get rid of two of the parameters.

I can Hessenberg reduce to get rid of (Subscript[γ, 13]) i.e. 
	(({{(Subscript[d, 1]), (Subscript[γ, 12]), (Subscript[γ, 13])}, {(Subscript[γ, 12]), (Subscript[d, 2]), (Subscript[γ, 23])}, {(Subscript[γ, 13]), (Subscript[γ, 23]), (Subscript[d, 3])}}))->(({{(Subscript[d, 1]), (Subscript[τ, 12]), 0}, {(Subscript[τ, 12]), (Subscript[d, 2]), (Subscript[τ, 23])}, {0, (Subscript[τ, 23]), (Subscript[d, 3])}}))
in which case my size restrictions are the same and I can do the same scaling. 
	(({{(Subscript[d, 1]), (Subscript[γ, 12]), (Subscript[γ, 13])}, {(Subscript[γ, 12]), (Subscript[d, 2]), (Subscript[γ, 23])}, {(Subscript[γ, 13]), (Subscript[γ, 23]), (Subscript[d, 3])}}))->(({{(Subscript[d, 1]), (Subscript[τ, 12]), 0}, {(Subscript[τ, 12]), (Subscript[d, 2]), (Subscript[τ, 23])}, {0, (Subscript[τ, 23]), 1-(Subscript[d, 1])-(Subscript[d, 2])}}))
with |(Subscript[τ, ij])|≤1/2  and 0≤(Subscript[d, i])≤1

vec={1,1,1};
vec.(({{(Subscript[d, 1]), (Subscript[τ, 12]), 0}, {(Subscript[τ, 12]), (Subscript[d, 2]), (Subscript[τ, 23])}, {0, (Subscript[τ, 23]), (Subscript[d, 3])}})).vec

Expand[{1,-1,-1}.B.{1,-1,-1}]

B={{d1,τ12,0},{τ12,d2,τ23},{0,τ23,1-d1-d2}};
MatrixForm[B]

J1.J2//MatrixForm

The product trig formulas
 cos(a)cos(b)=(1/2)(cos(a+b) + cos(a-b))
  sin(a)sin(b)=(1/2)(cos(a+b) - cos(a-b))
   sin(a)cos(b)=(1/2)(sin(a+b) + sin(a-b))
    cos(a)sin(b)=(1/2)(sin(a+b) - sin(a-b))
    
    with ϕp=θ1+θ2 ϕm=θ1-θ2  we have θ1=(1/2)(ϕ1+ϕ2) and θ2=(1/2)(ϕ1-ϕ2)
    
    (({{Cos[θ1], Cos[θ2] Sin[θ1], Sin[θ1] Sin[θ2]}, {-Sin[θ1], Cos[θ1] Cos[θ2], Cos[θ1] Sin[θ2]}, {0, -Sin[θ2], Cos[θ2]}}))->(({{Cos[(1/2)(ϕ2+ϕ2)], Cos[θ2] Sin[θ1], Sin[θ1] Sin[θ2]}, {-Sin[θ1], Cos[θ1] Cos[θ2], Cos[θ1] Sin[θ2]}, {0, -Sin[θ2], Cos[θ2]}}))

I can Jacobi out either of the off-diagonal entry.  The problem with this is it distinguishes one of the entries and makes an “if” in the logic. Probably better to Jacobi so that (Subscript[τ, 12])=(Subscript[τ, 23])

B={{d1,τ12,0},{τ12,d2,τ23},{0,τ23,d3}};
J1 = {
{Cos[θ1], Sin[θ1],0},
{-Sin[θ1], Cos[θ1],0},
{0,0,1}
};
J2= {
{1,0,0},
{0, Cos[θ2],Sin[θ2]},
{0,-Sin[θ2], Cos[θ2]}
};
Fred=Simplify[Transpose[J2].Transpose[J1].B.J1.J2];
MatrixForm[Fred]

Simplify[
Solve[
Simplify[
((Fred[[1,3]])/(Cos[θ2]))==0],
Tan[θ2]]
]

Simplify[Fred[[2,1]]-Fred[[3,2]]==0/.θ2->ArcTan[((2 τ23 Sin[θ1])/(2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]))]]

Together[%267]

Solve[Fred[[1,3]]==0,Sin[θ2]]

Simplify[Fred[[2,1]]-Fred[[3,2]]==0/.{
Sin[θ2]->((2 τ23 Cos[θ2] Sin[θ1])/(2 τ12 Cos[2 θ1]+d1 Sin[2 θ1]-d2 Sin[2 θ1]))
}]

Wilma=Simplify[
Numerator[
Together[(τ12 ((Cos[θ1])^2)+(d1-d2) Cos[θ1] Sin[θ1]-((2 τ23 Cos[θ2] (τ12 Cos[θ1]+d1 Sin[θ1]))/(2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]))+((2 τ23 Cos[θ2] Sin[θ1] (2 d3 τ12 Cos[2 θ1]+(d1 d3-d2 d3+((τ23)^2)) Sin[2 θ1]))/(((2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]))^2))+((Sin[θ1])^2) (-τ12+((2 ((τ23)^2))/(2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]))))]
]
]

Solve[Wilma==0,Sin[θ1]]

Solve[{
Fred[[1,3]]==0,
Fred[[2,1]]-Fred[[3,2]]==0
},
{θ1,θ2}]

θ2SolC=Simplify[Solve[-τ23 Cos[θ2] Sin[θ1]+(1/2) (2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]) Sin[θ2]==0,Sin[θ1]][[1]]]

Simplify[Solve[-τ23 Cos[θ2] Sin[θ1]+(1/2) (2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]) Sin[θ2]==0,θ1]]

θ2SolS=Simplify[Solve[-τ23 Cos[θ2] Sin[θ1]+(1/2) (2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]) Sin[θ2]==0,Sin[θ2]][[1]]]
θ2SolS=Simplify[Solve[-τ23 Cos[θ2] Sin[θ1]+(1/2) (2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]) Sin[θ2]==0,Cos[θ2]][[1]]]

Tan[θ2]=((2 τ23 Sin[θ1])/(2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]))

Barney=Simplify[Fred/.θ2->ArcTan[((2 τ23 Sin[θ1])/(2 τ12 Cos[2 θ1]+(d1-d2) Sin[2 θ1]))]]

FullSimplify[Barney[[1,2]],{Element[τ23,Reals]}]

f[λ_]=CharacteristicPolynomial[B,λ]
CoefficientList[f[λ],λ]//TableForm

Factor[d1 d2-(d1^2) d2-d1 (d2^2)-((τ12)^2)+d1 ((τ12)^2)+d2 ((τ12)^2)-d2 ((τ13)^2)+2 τ12 τ13 τ23-d1 ((τ23)^2)]

RootThing[γ0_,γ1_]=x/.Solve[
γ0+γ1 x + (x^2) - (x^3)==0,
x]//Simplify

ContourPlot[RootThing[γ0,γ1][[1]], 
{γ0,-2,2},
{γ1, -2,2}]
ContourPlot[RootThing[γ0,γ1][[2]], 
{γ0,-2,2},
{γ1, -2,2}]
ContourPlot[RootThing[γ0,γ1][[3]], 
{γ0,-2,2},
{γ1, -2,2}]

d1+ d2 + d3==

B=Randomreal

Eigenvectors

Simplify[Eigenvectors[B,1]/.Root[-d1 d2 d3+d3 ((τ12)^2)+d2 ((τ13)^2)-2 τ12 τ13 τ23+d1 ((τ23)^2)+(d1 d2+d1 d3+d2 d3-((τ12)^2)-((τ13)^2)-((τ23)^2)) #1+(-d1-d2-d3) ((#1)^2)+((#1)^3)&,1]->χ
][[1]]

A= {{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}};
AA=Flatten[A];
Ass=Table[Element[AA[[i]],Reals],{i,1, Length[AA]}];
{U,S,V}=SingularValueDecomposition[ A,1];
Simplify[S[[1,1]],Ass]

?ComplexExpand