Skip to content

Optimize matrix inverse in MGWR #73

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
HPDell wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Optimize matrix inverse in MGWR #73

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
390b39f
edit: add a benchmark
HPDell Oct 5, 2023
de8dbca
edit: test set bandwidth as initial value
HPDell Oct 5, 2023
e0ebbb9
edit: test add benchmark label
HPDell Oct 5, 2023
8c1e083
edit: optimise matrix inverse
HPDell Oct 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 11 additions & 11 deletions src/gwmodelpp/GWRMultiscale.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -426,17 +426,17 @@ vec GWRMultiscale::fitVarSerial(const vec &x, const vec &y, const uword var, mat
if (mHasHatMatrix)
{
mat ci, si;
S = mat(mHasHatMatrix ? nDp : 1, nDp, fill::zeros);
S = mat(nDp, nDp, fill::zeros);
for (uword i = 0; i < nDp ; i++)
{
GWM_LOG_STOP_BREAK(mStatus);
vec w = mSpatialWeights[var].weightVector(i);
mat xtw = trans(x % w);
mat xtwx = xtw * x;
mat xtwy = xtw * y;
double xtwx = as_scalar(xtw * x);
double xtwy = as_scalar(xtw * y);
try
{
mat xtwx_inv = inv_sympd(xtwx);
double xtwx_inv = 1.0 / xtwx;
betas.col(i) = xtwx_inv * xtwy;
ci = xtwx_inv * xtw;
si = x(i) * ci;
Expand Down Expand Up @@ -570,12 +570,12 @@ double GWRMultiscale::bandwidthSizeCriterionVarCVSerial(BandwidthWeight *bandwid
vec w = bandwidthWeight->weight(d);
w(i) = 0.0;
mat xtw = trans(mXi % w);
mat xtwx = xtw * mXi;
mat xtwy = xtw * mYi;
double xtwx = as_scalar(xtw * mXi);
double xtwy = as_scalar(xtw * mYi);
try
{
mat xtwx_inv = inv_sympd(xtwx);
vec beta = xtwx_inv * xtwy;
double xtwx_inv = 1.0 / xtwx;
double beta = xtwx_inv * xtwy;
double res = mYi(i) - det(mXi(i) * beta);
cv += res * res;
}
Expand Down Expand Up @@ -607,11 +607,11 @@ double GWRMultiscale::bandwidthSizeCriterionVarAICSerial(BandwidthWeight *bandwi
vec d = mSpatialWeights[var].distance()->distance(i);
vec w = bandwidthWeight->weight(d);
mat xtw = trans(mXi % w);
mat xtwx = xtw * mXi;
mat xtwy = xtw * mYi;
double xtwx = as_scalar(xtw * mXi);
double xtwy = as_scalar(xtw * mYi);
try
{
mat xtwx_inv = inv_sympd(xtwx);
double xtwx_inv = 1.0 / xtwx;
betas.col(i) = xtwx_inv * xtwy;
mat ci = xtwx_inv * xtw;
mat si = mXi(i) * ci;
Expand Down
1 change: 1 addition & 0 deletions test/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ set(SAMPLE_DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data)
add_definitions(-DSAMPLE_DATA_DIR="${SAMPLE_DATA_DIR}")

set(ADDON_SOURCES
data/londonhp.cpp
data/londonhp100.cpp
TerminateCheckTelegram.cpp
FileTelegram.cpp
Expand Down
45 changes: 45 additions & 0 deletions test/testGWRMultiscale.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "gwmodelpp/spatialweight/CRSDistance.h"
#include "gwmodelpp/spatialweight/BandwidthWeight.h"
#include "gwmodelpp/spatialweight/SpatialWeight.h"
#include "londonhp.h"
#include "londonhp100.h"
#include "TerminateCheckTelegram.h"
#include "FileTelegram.h"
Expand Down Expand Up @@ -336,3 +337,47 @@ TEST_CASE("Multiscale GWR: cancel")
}

}

TEST_CASE("Large Data Repeate", "[benchmark]")
{
mat londonhp_coord, londonhp_data;
vector<string> londonhp_fields;
if (!read_londonhp(londonhp_coord, londonhp_data, londonhp_fields))
{
FAIL("Cannot load londonhp100 data.");
}

vec y = londonhp_data.col(0);
mat x = join_rows(ones(londonhp_data.n_rows), londonhp_data.cols(uvec({1, 3})));
uword nVar = 3;

BENCHMARK("specified bandwidth")
{
vector<SpatialWeight> spatials;
vector<bool> preditorCentered;
vector<GWRMultiscale::BandwidthInitilizeType> bandwidthInitialize;
vector<GWRMultiscale::BandwidthSelectionCriterionType> bandwidthSelectionApproach;
for (size_t i = 0; i < nVar; i++)
{
CRSDistance distance;
BandwidthWeight bandwidth(64, true, BandwidthWeight::Bisquare);
spatials.push_back(SpatialWeight(&bandwidth, &distance));
preditorCentered.push_back(i != 0);
bandwidthInitialize.push_back(GWRMultiscale::BandwidthInitilizeType::Initial);
bandwidthSelectionApproach.push_back(GWRMultiscale::BandwidthSelectionCriterionType::CV);
}

GWRMultiscale algorithm;
algorithm.setCoords(londonhp_coord);
algorithm.setDependentVariable(y);
algorithm.setIndependentVariables(x);
algorithm.setSpatialWeights(spatials);
algorithm.setHasHatMatrix(true);
algorithm.setPreditorCentered(preditorCentered);
algorithm.setBandwidthInitilize(bandwidthInitialize);
algorithm.setBandwidthSelectionApproach(bandwidthSelectionApproach);
algorithm.setBandwidthSelectThreshold(vector(3, 1e-5));
algorithm.fit();
return algorithm.betas();
};
}