Simple J4 Implementation

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}
]
]
HamiltonJacobi4Sym[S_]:=Module[
{T=ConstantArray[0,{3,3}],u,σ,v},
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]]);
{u,σ,v}=SingularValueDecomposition[T,1];{u,σ,v}=Map[Flatten,{u,σ,v}];
u=Prepend[u,0];
v=Prepend[v,0];
u[[2]]+=1;
v[[2]]+=1;
u=u/Norm[u];
v=v/Norm[v];
HJtensorProduct[u,v]
]

Clear[HamiltonJacobi4Skw]
HamiltonJacobi4Skw[K_]:=Module[
{p,q,i=Quaternion[0,1,0,0],pQ,qQ,Np,Nq},
(* The following are Niloufer M 6.6 and 6.7 *)
p=0.5 {-(K[[1,2]]+K[[3,4]]),(K[[2,4]]-K[[1,3]]),-(K[[1,4]]+K[[2,3]])};
q= 0.5{   K[[1,2]]-K[[3,4]],     K[[2,4]]+K[[1,3]],      K[[1,4]]-K[[2,3]]};
{pQ,qQ}={Apply[Quaternion,Join[{0},p]],Apply[Quaternion,Join[{0},q]]};
(* The stuff below does not work *)
(* Not sure why but it almost has to be in the Tensor Product*)
(*{p,q}={Apply[Quaternion,p],Apply[Quaternion,q]};
{x,y}={((Abs[p]- i**p)/(Abs[Abs[p]-i**p])),((Abs[q]- i **q)/(Abs[Abs[q]-i **q]))};
{x,y}={Apply[List,x],Apply[List,y]};
HJtensorProduct[x,y]*)
(* direct implementation *)
{Np,Nq}=Map[Norm,{p,q}];
Qp=(1/(Abs[Abs[pQ]-i**pQ]))(({{Np+p[[1]], 0, -p[[3]], p[[2]]}, {0, Np+p[[1]], p[[2]], p[[3]]}, {p[[3]], -p[[2]], Np+p[[1]], 0}, {-p[[2]], -p[[3]], 0, Np+p[[1]]}}));
Qq=(1/(Abs[Abs[qQ]-i**qQ]))(({{Nq+q[[1]], 0, q[[3]], -q[[2]]}, {0, Nq+q[[1]], q[[2]], q[[3]]}, {-q[[3]], -q[[2]], Nq+q[[1]], 0}, {q[[2]], -q[[3]], 0, Nq+q[[1]]}}));
Qp.Qq
]

NonUpperTriangularity[A_List]:=Module[
{AUpper},
AUpper=UpperTriangularize[ A];
1-Norm[AUpper,"Frobenius"]/Norm[A,"Frobenius"]
]
NonUpperTriangularity[A_SparseArray]:=Module[
{AUpper},
AUpper=UpperTriangularize[ A];
1-Norm[AUpper,"Frobenius"]/Norm[A,"Frobenius"]
]

Testing

n=4;
AOriginal=RandomReal[{-1,1},{n,n}];
ASym= AOriginal+Transpose[AOriginal];
S=ASym;
HSym=HamiltonJacobi4Sym[S];
S=Chop[HSym.S.Transpose[HSym]];
MatrixPlot[S,
PlotLabel->Norm[HSym.Transpose[HSym]-IdentityMatrix[4]]
]

n=4;
AOriginal=RandomReal[{-1,1},{n,n}];
ASkw= AOriginal-Transpose[AOriginal];
K=ASkw;
HSkw=HamiltonJacobi4Skw[K];
K=Chop[HSkw.K.Transpose[HSkw]];
MatrixPlot[K,
PlotLabel->Norm[HSkw.Transpose[HSkw]-IdentityMatrix[4]]
]

Together

n=4;
AOriginal=RandomReal[{-1,1},{n,n}];
A=AOriginal;
{S,K}=0.5{ A+Transpose[A],A-Transpose[A]};
HSym=HamiltonJacobi4Sym[S];
HSkw=HamiltonJacobi4Skw[K];
{MatrixPlot[Chop[HSym.S.Transpose[HSym]],
PlotLabel->Norm[HSym.Transpose[HSym]-IdentityMatrix[4]]
],
MatrixPlot[Chop[HSkw.K.Transpose[HSkw]],
PlotLabel->Norm[HSkw.Transpose[HSkw]-IdentityMatrix[4]]
]
}

n=4;
AOriginal=RandomReal[{-1,1},{n,n}];
A=AOriginal;
{S,K}=0.5{ A+Transpose[A],A-Transpose[A]};
HSym=HamiltonJacobi4Sym[S];
S=Chop[HSym.S.Transpose[HSym]];
KTemp=Chop[HSym.K.Transpose[HSym]]
HSkw=HamiltonJacobi4Skw[KTemp];
HBig=HSkw.HSym;
{MatrixPlot[Chop[HSym.S.Transpose[HSym]],
PlotLabel->Norm[HSym.Transpose[HSym]-IdentityMatrix[4]]
],
MatrixPlot[Chop[HBig.K.Transpose[HBig]],
PlotLabel->Norm[HBig.Transpose[HBig]-IdentityMatrix[4]]
],
MatrixPlot[ANew=Chop[HBig.A.Transpose[HBig]],
PlotLabel->Norm[HBig.Transpose[HBig]-IdentityMatrix[4]]
]
}

MatrixForm[ANew]

Quaternion Basis

The Quaternion basis is the secret for reinterpreting an objective function for the matrices in 
terms of a pair of quaternions.This shouldbe pretty simple.  First lets do the Frobenius norm 
of the matrix.

QBasis=(({{(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}})), (({{0, 1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0}})), (({{0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0, 0}})), (({{0, 0, 0, 1}, {0, 0, -1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}}))}, {(({{0, -1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, -1}, {0, 0, 1, 0}})), (({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}})), (({{0, 0, 0, -1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {-1, 0, 0, 0}})), (({{0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}, {0, 1, 0, 0}}))}, {(({{0, 0, -1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}, {0, -1, 0, 0}})), (({{0, 0, 0, 1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {1, 0, 0, 0}})), (({{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, -1}})), (({{0, -1, 0, 0}, {-1, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}}))}, {(({{0, 0, 0, -1}, {0, 0, -1, 0}, {0, 1, 0, 0}, {1, 0, 0, 0}})), (({{0, 0, -1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, 1, 0, 0}})), (({{0, 1, 0, 0}, {1, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}})), (({{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, 1}}))}}));
αMat=Sum[α[i,j]QBasis[[i,j]],{j,4},{i,4}];
Perm={13,9,14,5,10,15,1,6,11,16,2,7,12,3,8,4};
Eqs=Flatten[
Table[
a[i,j]==αMat[[i,j]],{j,4},{i,4}]
]
Vars=Flatten[Array[α,{4,4}]]
PVars=Permute[Vars,Perm]

A=Flatten[Sum[ α[i,j] QBasis[[i,j]],{i,4},{j,4}]];
Vars=Flatten[Table[α[i,j],{i,4},{j,4}]];
PVars=Permute[Vars,Perm];
MatrixForm[M = (1/2)CoefficientArrays[A,Vars][[2]]];
MatrixForm[PM =Permute[ (1/2)CoefficientArrays[A,PVars][[2]],Perm]];
(* Note M is orthogonal in fact a rotation and involutary *)
(* Det[M]=1 and M.M=Id *)
(* ½(M+I) is idempotent *)
Map[MatrixForm,{M.M,PM.PM}]
Map[MatrixPlot,{M,PM}]

MatrixPlot[PM]

So the connection between the As and the αs is given by 2M

MatrixForm[AToα=Simplify[Partition[ 2 M.RHS,4]]]
MatrixForm[αToA=Simplify[Partition[ 2 M.Vars,4]]]

I am actually interested in setting the most off-diagonal parts below the diagonal to zero.

MatrixForm[
MSub=CoefficientArrays[{
αToA[[1,4]]==0,
αToA[[1,3]]==0,
αToA[[2,4]]==0,
αToA[[1,2]]==0,
αToA[[2,3]]==0
},
Vars
][[2]]
]

Drop[M,{13,}]

PM=Permute[M,{13,9,14,5,10,15,1,6,11,16,2,7,12,3,8,4}]
MatrixForm[PM.PM]

MatrixPlot[M]

So it is easy to look for

{MatrixForm[PseudoInverse[1/2 MSub]],MatrixForm[1/2 Transpose[MSub]]}

FoldM=Partition[M,{4,4}];
OneWay=Simplify[FoldM.Array[a,{4,4}]];
Dimensions[OneWay]

Dimensions[M]

MatrixForm[
NSpace=NullSpace[MSub]
]
Dimensions[NSpace]

TableForm[{
{
αToA[[1,4]],
αToA[[1,3]],
αToA[[2,4]]
},
Normal[2M[[1]]]
}
]

CoefficientLi{
αToA[[1,4]],
αToA[[1,3]],
αToA[[2,4]]
}

Looking at the matrices reveals something interesting. With this partioning m1 = 4 Invm1. Why does this happen?

m1-4 Invm1

Why does this happen?

Diagonal Entries in 4×4s

The combinations that give clean diagonal output are pretty easy to find and give the fours.

QDD={QBasis[[1,1]],QBasis[[2,2]],QBasis[[3,3]],QBasis[[4,4]]};
ADD=1/2{
{1,1,1,1},
{1,1,-1,-1},
{1,-1,1,-1},
{1,-1,-1,1}
};
Map[MatrixForm,ADD.QDD]
MatrixForm[ADD.ADD]

Off-Diagonal Entries in Diagonal 2×2s inside the 4×4s

The combinations that gives the off diagonal entries in the diagonal 2×2s are pretty easy to find and also give fours.

QOD={QBasis[[1,2]],QBasis[[2,1]],QBasis[[3,4]],QBasis[[4,3]]};
AOD=1/2{
{1,-1,-1,1},
{-1,1,-1,1},
{1,1,1,1},
{-1,-1,1,1}
};

Map[MatrixForm,AOD.QOD]
MatrixForm[ AOD.AOD]

Diagonal Entries in Off-Diagonal 2×2s inside the 4×4s

The combinations that gives the diagonal entries in the off-diagonal 2×2s are pretty easy to find and also give fours.

QDO={QBasis[[1,3]],QBasis[[2,4]],QBasis[[3,1]],QBasis[[4,2]]};
ADO=1/2{
{-1,1,1,-1},
{-1,1,-1,1},
{1,1,1,1},
{1,1,-1,-1}
};
Map[MatrixForm,ADO.QDO]
MatrixForm[ADO.ADO]

Off-Diagonal Entries in Off-Diagonal 2×2s inside the 4×4s

The combinations that gives the off-diagonal entries in the off-diagonal 2×2s are pretty easy to find and also give fours.

QBOO={QBasis[[1,4]],QBasis[[2,3]],QBasis[[3,2]],QBasis[[4,1]]};
AOO=1/2{
{-1,1,1,-1},
{1,-1,1,-1},
{1,1,1,1},
{-1,-1,1,1}
};
Map[MatrixForm,AOO.QBOO
]
MatrixForm[AOO.AOO]

Optimization Again

I have a function (can be quadratic) of two 4 dimensional vector variables.  It amost certainly has some symmetrires. I have no idea what they might be?

f[p_,q_]:= p.A.p

p=Quaternion[p1,p2,p3,p4];
x=Quaternion[x1,x2,x3,x4];
conjThing=Simplify[x**p**Conjugate[x]]/.(x1^2)+(x2^2)+(x3^2)+(x4^2)->1;
conjThing=Apply[List,conjThing];
MatrixForm[Simplify[ D[conjThing,{{x1,x2,x3,x4},1}]]]
MatrixForm[conjThing]

NewConjThing=Simplify[x**p**Conjugate[x]];
NewConjThing=Apply[List,NewConjThing];
MatrixForm[NewConjThing]

May 14th

Here is the question.  What happens if I do the diagonalization trick on the symmetric part of the matrix.  I am going to deconstruct the original J4 code and experiment.

Code talk through:
The first portion translates the (assumed input symmetric matrix S) into the T matrix. 

T=1/4(({{S[[1,1]]+S[[2,2]]-S[[3,3]]-S[[4,4]], 2(S[[2,3]]-S[[1,4]]), 2(S[[2,4]]+S[[1,3]])}, {2(S[[2,3]]+S[[1,4]]), S[[1,1]]-S[[2,2]]+S[[3,3]]-S[[4,4]], 2(S[[3,4]]-S[[1,2]])}, {2(S[[2,4]]-S[[1,3]]), 2(S[[3,4]]+S[[1,2]]), S[[1,1]]-S[[2,2]]-S[[3,3]]+S[[4,4]]}}))
when I fix this I need to undo the symmetrizing that has occured.  In other words half of the S[[2,3]]s are S[[3,2]] etc when the formula is desymmetrized. The easiest
way to do this is to use the linearity of the representation

M4ToQxQ[A_]:= Module[
{S,K,cS,pS,qS,rS,,u,σ,v,pK,qK},
S=0.5(A + Transpose[A]);
K=A-A;
cS=0.25 Trace[S];
{pS,qS,rS}=0.25(({{S[[1,1]]+S[[2,2]]-S[[3,3]]-S[[4,4]], 2(S[[2,3]]-S[[1,4]]), 2(S[[2,4]]+S[[1,3]])}, {2(S[[2,3]]+S[[1,4]]), S[[1,1]]-S[[2,2]]+S[[3,3]]-S[[4,4]], 2(S[[3,4]]-S[[1,2]])}, {2(S[[2,4]]-S[[1,3]]), 2(S[[3,4]]+S[[1,2]]), S[[1,1]]-S[[2,2]]-S[[3,3]]+S[[4,4]]}}));
(* Modified to rerurn all three SVs *)
{u,σ,v}=SingularValueDecomposition[{pS,qS,rS},3];
Print[MatrixForm[Chop[Transpose[u].{pS,qS,rS}.v]]];
{pK,qK}=0.5(({{-(K[[1,2]]+K[[3,4]]), K[[2,4]]-K[[1,3]], -(K[[1,4]]+K[[2,3]])}, {  (K[[1,2]]-K[[3,4]]), K[[2,4]]+K[[1,3]],    K[[1,4]]-K[[2,3]]}}));

]

n=4;
A=RandomReal[{-1,1},{n,n}];
M4ToQxQ[A]

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}
]
]

HamiltonJacobi4Sym[S_]:=Module[
{T=ConstantArray[0,{3,3}],u,σ,v},
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;
{u,σ,v}=SingularValueDecomposition[T,1];
{u,σ,v}=Map[Flatten,{u,σ,v}];
u=Prepend[u,0];
v=Prepend[v,0];
u[[2]]+=1;
v[[2]]+=1;
u=u/Norm[u];
v=v/Norm[v];
HJtensorProduct[u,v]
]

n=4;
AOriginal=RandomReal[{-1,1},{n,n}];
ASym= AOriginal+Transpose[AOriginal];
S=ASym;
M4toQxQ[S]

HSym=HamiltonJacobi4Sym[S];
S=Chop[HSym.S.Transpose[HSym]];
MatrixPlot[S,
PlotLabel->Norm[HSym.Transpose[HSym]-IdentityMatrix[4]]
]

MatrixForm[TOut]

Optimization

I am trying to optimize a quadratic in three variables subject to two very clean quadratic equality 
constraints.  So

A= RandomReal[{-1,1},{8,8}];
A=A.Transpose[A];
A=A+Transpose[A];
α=Join[Array[α1,4],Array[α2,4]];
B1=SparseArray[Band[{1,1},{4,4}]->1,{8,8}];
B2=SparseArray[Band[{5,5}]->1,{8,8}]
NMinimize[{α.A.α,{α.B1.α==1,α.B2.α==1}},α]

A=RandomReal[{-1,1},{8,8}];
A=A.Transpose[A];
A=A+Transpose[A];
α=Array[α1,{8}];
{val,Subs}=Minimize[{ (1/2)α.A.α, Norm[Take[α,4]]==1, Norm[Drop[α,4]]==1},α]

{val-(1/2)α.A.α, Norm[Take[α,4]]-1, Norm[Drop[α,4]]-1}/.Subs

Break

Clear[n,A,AOriginal,Id,H,HBig]
n=12;
AOriginal= SparseArray[{
Band[{1,1}]->RandomReal[{-1,1},n],
Band[{1,2}]->RandomReal[{-1,1},n-1],
Band[{1,3}]->RandomReal[{-1,1},n-2],
Band[{1,3}]->RandomReal[{-1,1},n-3]
},
{n,n}];
AOriginal=AOriginal+Transpose[AOriginal];
A= AOriginal;
Id=IdentityMatrix[n];

Join[
{MatrixPlot[A]},
Table[
ASub=A[[i;;i+3,i;;i+3]];
H=HamiltonJacobi4[ASub];
HBig=Id;
HBig[[i;;i+3,i;;i+3]]=H;
A=Chop[HBig.A.Transpose[HBig]];
MatrixPlot[A,PlotLabel->Map[NonUpperTriangularity,{A,AOriginal}]],
{i,1,n-3}]
]
Transpose[Map[Eigenvalues,{A,AOriginal}]]

ListPlot[
Map[Abs[Diagonal[#]]&,{A,AOriginal}],
Joined->True,
PlotStyle->{Red,Blue},
Mesh->All]

Norm[Eigenvalues[A]-Eigenvalues[AOriginal]]

Understanding the pattern

ASym J4 Implementation

The decomposition
	A=(Subscript[A, sym])+(Subscript[A, skw])
always works.

n=4;
AOriginal=RandomReal[{-1,1},{n,n}];
A=AOriginal;
{ASym,ASkw}={(A+Transpose[A])/2,(A-Transpose[A])/2};
H=HamiltonJacobi4[ASym];
ASkwNew=Chop[H.ASkw.Transpose[H]];
MatrixPlot[ASkwNew]
MatrixPlot[Chop[H.ASym.Transpose[H]]]

Testing

n=4;
AOriginal=RandomReal[{-1,1},{n,n}];
AOriginal= AOriginal+Transpose[AOriginal];
A=AOriginal;
H=HamiltonJacobi4[A];
A=Chop[H.A.Transpose[H]];
Transpose[Map[Eigenvalues,{A,AOriginal}]]
MatrixPlot[A,
PlotLabel->Map[NonUpperTriangularity,{A,AOriginal}]]

ListPlot[Map[Abs[Diagonal[#]]&,{A,AOriginal}],
Joined->True,
PlotStyle->{Red,Blue},
Mesh->All]

testing on a band matrix

Clear[n,A,AOriginal,Id,H,HBig]
n=12;
AOriginal= SparseArray[{
Band[{1,1}]->RandomReal[{-1,1},n],
Band[{1,2}]->RandomReal[{-1,1},n-1],
Band[{1,3}]->RandomReal[{-1,1},n-2],
Band[{1,3}]->RandomReal[{-1,1},n-3]
},
{n,n}];
AOriginal=AOriginal+Transpose[AOriginal];
A= AOriginal;
Id=IdentityMatrix[n];

Join[
{MatrixPlot[A]},
Table[
ASub=A[[i;;i+3,i;;i+3]];
H=HamiltonJacobi4[ASub];
HBig=Id;
HBig[[i;;i+3,i;;i+3]]=H;
A=Chop[HBig.A.Transpose[HBig]];
MatrixPlot[A,PlotLabel->Map[NonUpperTriangularity,{A,AOriginal}]],
{i,1,n-3}]
]
Transpose[Map[Eigenvalues,{A,AOriginal}]]

ListPlot[
Map[Abs[Diagonal[#]]&,{A,AOriginal}],
Joined->True,
PlotStyle->{Red,Blue},
Mesh->All]

Norm[Eigenvalues[A]-Eigenvalues[AOriginal]]

Understanding the pattern

Clear[n,A,AOriginal,Id,H,HBig]
n=12;
AOriginal= RandomReal[{-1,1},{n,n}];
AOriginal=AOriginal.Transpose[AOriginal];
A= AOriginal;
Id=IdentityMatrix[n];

range={1,2,3,4};
ASub=A[[range,range]];
H=HamiltonJacobi4[ASub];
HBig=Id;
HBig[[range,range]]=H;
A=Chop[HBig.A.Transpose[HBig]];
MatrixPlot[A,PlotLabel->NonUpperTriangularity[A]]
MatrixPlot[HBig]

In H.A phase uses the pivot rows

MatrixForm[H]

Again

Clear[n,A,AOriginal,Id,H,HBig]
n=12;
AOriginal= RandomReal[{-1,1},{n,n}];
AOriginal=AOriginal.Transpose[AOriginal];
A= AOriginal;
Id=IdentityMatrix[n];

Table[
range={1,2,i,i+1};
ASub=A[[range,range]];
H=HamiltonJacobi4[ASub];
HBig=Id;
HBig[[range,range]]=H;
A=Chop[HBig.A.Transpose[HBig]];
MatrixPlot[A,PlotLabel->NonUpperTriangularity[A]],
{i,3,11,2}]

J8 Start

I should be able to assemble multiple J4s into a suitable J8

SetAttributes[HJ,HoldFirst]
 HJ[A_,range_]:= Module[{ASub=A[[range,range]]},
Switch[Length[range],
4, HamiltonJacobi4[ASub],
8, HamiltonJacobi8[ASub]
]
]

n=8;
AOriginal= RandomReal[{-1,1},{n,n}];
AOriginal=AOriginal+Transpose[AOriginal];
A= AOriginal;

HTemp=ConstantArray[0,{n,n}];
ranges={{1,2,3,4},{5,6,7,8}};
Do[
r=ranges[[i]];
HTemp[[r,r]]=HamiltonJacobi4[A[[r,r]]],
{i,1,2}]
A=Chop[HTemp.A.Transpose[HTemp]];
MatrixPlot[A,PlotLabel->NonUpperTriangularity[A]]

HTemp=ConstantArray[0,{n,n}];
ranges={{1,2,5,6},{3,4,7,8}};
Do[
r=ranges[[i]];
HTemp[[r,r]]=HamiltonJacobi4[A[[r,r]]],
{i,1,2}]
A=Chop[HTemp.A.Transpose[HTemp]];
MatrixPlot[A,PlotLabel->NonUpperTriangularity[A]]

HTemp=ConstantArray[0,{n,n}];
ranges={{1,2,7,8},{3,4,5,6}};
Do[
r=ranges[[i]];
HTemp[[r,r]]=HamiltonJacobi4[A[[r,r]]],
{i,1,2}]
A=Chop[HTemp.A.Transpose[HTemp]];
MatrixPlot[A,PlotLabel->NonUpperTriangularity[A]]

MatrixForm[A]

MmaEvals=Eigenvalues[AOriginal];
((Eigenvalues[A]-MmaEvals)/MmaEvals)

A= AOriginal;
Id=IdentityMatrix[n];

range={1,2,5,6};
ASub=A[[range,range]];
H=HamiltonJacobi[ASub];
HBig=Id;
HBig[[range,range]]=H;
A=Chop[HBig.A.Transpose[HBig]];
MatrixPlot[A,PlotLabel->NonUpperTriangularity[A]]
ListLogPlot[Abs[((Eigenvalues[A]-Eigenvalues[AOriginal])/(Eigenvalues[AOriginal]))],
PlotRange->All
]

ListPlot[
{Abs[Diagonal[A]],Abs[Eigenvalues[A]]},
Joined->True,
PlotStyle->{Red,Green},
Mesh->All]

ListPlot[Map[Abs[Diagonal[#]]&,{A,AOriginal}],
Joined->True,
PlotStyle->{Red,Blue},
Mesh->All]