diff --git a/docs/site/builtins-reference.md b/docs/site/builtins-reference.md index 22b335866cb..e99182a87ca 100644 --- a/docs/site/builtins-reference.md +++ b/docs/site/builtins-reference.md @@ -28,6 +28,7 @@ limitations under the License. * [`tensor`-Function](#tensor-function) * [DML-Bodied Built-In functions](#dml-bodied-built-in-functions) * [`confusionMatrix`-Function](#confusionmatrix-function) + * [`coresetDT`-Function](#coresetdt-function) * [`correctTypos`-Function](#correcttypos-function) * [`cspline`-Function](#cspline-function) * [`csplineCG`-Function](#csplineCG-function) @@ -207,6 +208,63 @@ y = toOneHot(X, numClasses) [ConfusionSum, ConfusionAvg] = confusionMatrix(P=z, Y=y) ``` +## `coresetDT`-Function + +Coreset selection for Classification via a depth-bounded decision tree. + +The method is adapted from the datamap-driven coreset approach of Hadar et al., +"Datamap-Driven Tabular Coreset Selection for Classifier Training" +(https://www.vldb.org/pvldb/vol18/p876-razmadze.pdf). + +A decision tree is trained on (X, y) and each leaf defines a region which +is classified by its majority fraction and size. + +majority_fraction < psi -> hard (mixed labels, keep whole)
+majority_fraction >= psi, n <= tau -> ambiguous (pure but small, keep whole)
+majority_fraction >= psi, n > tau -> easy (large + pure, sample down) + +Hard and ambiguous regions are kept in full. Easy regions are sampled down +to samp_ratio, processed largest-first, until the target size is reached, +keeping any unvisited regions whole once they fit. The result cannot +shrink below:
+ n_min ~ |hard rows| + |ambiguous rows| + samp_ratio * |easy rows|
+and requests below this floor return roughly the floor (with a warning). + +### Usage + +```r +[Xc, yc] = coresetDT(X, y, fraction, samp_ratio, psi, tau, max_depth, min_leaf, seed, verbose) +``` + +### Arguments + +| Name | Type | Default | Description | +| :--------- | :------------- | :------ | :---------- | +| X | Matrix[Double] | --- | Feature matrix in recoded/binned representation (as required by decisionTree)| +| y | Matrix[Double] | --- | Continuous label vector (as required by decisionTree)| +| fraction | Double | 0.3 | Target coreset size | +| samp_ratio | Double | 0.03 | Thinning ratio applied to easy regions | +| psi | Double | 0.90 | Homogeneity threshold separating hard regions | +| tau | Int | 10 | Absolute size threshold separating easy from ambiguous regions | +| max_depth | Int | 12 | Tree depth for the datamap | +| min_leaf | Int | 5 | Minimum samples per leaf | +| seed | Int | -1 | Seed for thinning| +| verbose | Boolean | FALSE | Print information and summary | + +### Returns + +| Name | Type | Description | +| :--- | :------------- | :---------- | +| Xc | Matrix[Double] | Coreset feature matrix | +| yc | Matrix[Double] | Coreset label vector | + +### Example + +```r +X = read($X) +y = read($y) +[Xc, yc] = coresetDT(X = X, y = y, fraction = 0.1, verbose = TRUE) +``` ## `correctTypos`-Function diff --git a/scripts/builtin/coresetDT.dml b/scripts/builtin/coresetDT.dml new file mode 100644 index 00000000000..c6f7b37ebf6 --- /dev/null +++ b/scripts/builtin/coresetDT.dml @@ -0,0 +1,198 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +# Coreset selection for Classification via a depth-bounded decision tree. +# +# The method is adapted from the datamap-driven coreset approach of Hadar et al., +# "Datamap-Driven Tabular Coreset Selection for Classifier Training" +# (https://www.vldb.org/pvldb/vol18/p876-razmadze.pdf). +# +# A decision tree is trained on (X, y) and each leaf defines a region which +# is classified by its majority fraction and size. +# +# majority_fraction < psi -> hard (mixed labels, keep whole) +# majority_fraction >= psi, n <= tau -> ambiguous (pure but small, keep whole) +# majority_fraction >= psi, n > tau -> easy (large + pure, sample down) +# +# Hard and ambiguous regions are kept in full. Easy regions are sampled down +# to samp_ratio, processed largest-first, until the target size is reached, +# keeping any unseen regions whole once they fit. The result cannot +# shrink below: +# n_min ~ |hard rows| + |ambiguous rows| + (samp_ratio * |easy rows|) +# and requests below this floor return roughly the floor (with a warning). +# +# INPUT: +# ------------------------------------------------------------------------------ +# X Feature matrix in recoded/binned representation +# (as required by decisionTree) +# y Continuous label vector (as required by decisionTree) +# fraction Target coreset size as a fraction of nrow(X) (default 0.3) +# samp_ratio Thinning ratio applied to easy regions (default 0.03) +# psi Homogeneity threshold separating hard regions (default 0.90) +# tau Absolute size threshold separating easy from ambiguous +# regions among pure leaves (default 10) +# min_leaf Minimum samples per leaf (default 5) +# max_depth Bounded tree depth for the datamap (default 12) +# seed Seed for reproducible thinning, -1 for random (default -1) +# verbose Print information and summary (default FALSE) +# ------------------------------------------------------------------------------ +# +# OUTPUT: +# ------------------------------------------------------------------------------ +# Xc Coreset feature matrix +# yc Coreset label vector +# ------------------------------------------------------------------------------ + +m_coresetDT = function(Matrix[Double] X, Matrix[Double] y, + Double fraction = 0.3, Double samp_ratio = 0.03, Double psi = 0.90, Int tau = 10, + Int max_depth = 12, Int min_leaf = 5, + Int seed = -1, Boolean verbose = FALSE) + + return(Matrix[Double] Xc, Matrix[Double] yc) +{ + n = nrow(X); + k = as.integer(max(y)); + nTarget = round(fraction * n) + + # Step 1: build the datamap + ctypes = cbind(matrix(1, 1, ncol(X)), matrix(2, 1, 1)); + tree = decisionTree(X=X, y=y, ctypes=ctypes, max_depth=max_depth, + min_leaf=min_leaf, min_split=2*min_leaf, max_features=1.0, seed=seed); + + # route every row to its leaf (tree traversal adapted from decisionTreePredict.dml) + numNodes = as.integer(ncol(tree) / 2); + nodes = matrix(tree, rows=numNodes, cols=2); + features = nodes[,1]; # split feature per node, 0 for leaves + values = nodes[,2]; # split value per node, class label for leaves + + Tidx = matrix(1, n, 1); + noChange = FALSE; + i = 1; + while(!noChange & i <= max_depth) { + P = table(seq(1,n), Tidx, n, numNodes); # one hot row->node + f = P %*% features; # split feature per row + v = P %*% values ; # split value per row + isNotLeaf = f > 0; # check if internal node + xv = rowSums(X * table(seq(1,n), max(f,1), n, ncol(X))); # xv[r] = X[r, f[r]] rows value for current split feature + Tidx_new = ifelse(isNotLeaf, 2*Tidx + (xv > v), Tidx); # descend to child + noChange = (sum(Tidx != Tidx_new) == 0); + Tidx = Tidx_new; + i = i + 1; + } + P = table(seq(1,n), Tidx, n, numNodes); + + # region stats + C = table(Tidx, y, numNodes, k); # label counts per region + s = rowSums(C); # region sizes + mf = rowMaxs(C) / max(s, 1); # majority fraction + exists = (s > 0); + isHard = exists & (mf < psi); + isAmb = exists & (mf >= psi) & (s <= tau); + isEasy = exists & (mf >= psi) & (s > tau); + + nHardRows = sum(isHard * s); + nAmbRows = sum(isAmb * s); + nEasyRows = sum(isEasy * s); + nFloor = nHardRows + nAmbRows + ceil(samp_ratio * nEasyRows); + + if(verbose) { + pHard = round(1000 * nHardRows / n) / 10; + pAmb = round(1000 * nAmbRows/ n) / 10; + pEasy = round(1000 * nEasyRows / n) / 10; + + print("coresetDT: datamap (psi=" + psi + ", tau=" + tau + ", max_depth=" + max_depth + + ", samp_ratio=" + samp_ratio + ") over " + n + " rows"); + + print("coresetDT: hard (mixed, mf<" + psi + "): " + sum(isHard) + " regions, " + + nHardRows + " rows (" + pHard + "% of data) -> kept whole"); + print("coresetDT: ambiguous (pure, size<=" + tau + "): " + sum(isAmb) + " regions, " + + nAmbRows + " rows (" + pAmb + "% of data) -> kept whole"); + print("coresetDT: easy (pure, size>" + tau + "): " + sum(isEasy) + " regions, " + + nEasyRows + " rows (" + pEasy + "% of data) -> thinned to "+ samp_ratio); + + print("coresetDT: size floor n_min = " + nFloor + " (" + (round(1000 * nFloor/ n) / 10) + "% of n)"); + } + + if(nTarget < nFloor) + print("coresetDT: WARN: requested size " + nTarget + " is below the floor " + nFloor + ". Returning ~ " + nFloor + + " rows. Lower samp_ratio, raise psi, or loosen the tree to go smaller."); + + # Step 2: select the coreset + keep = P %*% (isHard + isAmb); + cur = sum(keep); + + # Spend remaining budget on easy regions + nEasy = as.integer(sum(isEasy)); + + if(cur < nTarget) { + easyIdx = removeEmpty(target=seq(1,numNodes), margin="rows", select=isEasy); + easySize = removeEmpty(target=s, margin="rows", select=isEasy); + + ordered = order(target=easySize, by=1, decreasing=TRUE, index.return=TRUE); + randVals = rand(rows=n, cols=1, seed=seed); + visited = matrix(0, numNodes, 1); + remaining = nEasyRows; + + done = FALSE; + for(i in 1:nEasy) { + if(!done) { + pos = as.scalar(ordered[i,1]); + regid = as.scalar(easyIdx[pos,1]); + regSize = as.scalar(easySize[pos,1]); + + visited[regid,1] = 1; + remaining = remaining - regSize; + + # thin region to samp_ratio, its the k smallest random scores within region + keeprows = max(round(samp_ratio * regSize), 1); + I = P[,regid]; # rows in this region + regRows = removeEmpty(target=seq(1,n), margin="rows", select=I); # global ids of region rows + regScores = removeEmpty(target=randVals, margin="rows", select=I); # their random scores + rank = order(target=regScores, by=1, index.return=TRUE); # rank rows within region + + chosenIds = table(seq(1,keeprows), rank[1:keeprows,1], keeprows, nrow(regRows)) %*% regRows; + keep = keep + table(chosenIds, 1, n, 1); + cur = cur + keeprows; + + # if all unvisited easy regions fit into the budget, add them + if(cur + remaining <= nTarget) { + keep = keep + P %*% (isEasy *(1 - visited)); + cur = cur + remaining; + done = TRUE; + } + } + } + } + + # extract coreset rows + toKeepIdx = removeEmpty(target=seq(1,n), margin="rows", select=keep); + Psel = table(seq(1,nrow(toKeepIdx)), toKeepIdx, nrow(toKeepIdx), n); + Xc = Psel %*% X; + yc = Psel %*% y; + + if(verbose) { + sel = nrow(Xc); + mustKeep = nHardRows + nAmbRows; + sampledEasy = sel - mustKeep; + print("coresetDT: selected " + sel + " / " + n + " rows (" + (round(1000 * sel / n) / 10) + "% of data) " + mustKeep + + " must-keep (hard+ambiguous) + " + sampledEasy + " sampled from easy"); + } +} diff --git a/src/main/java/org/apache/sysds/common/Builtins.java b/src/main/java/org/apache/sysds/common/Builtins.java index f5719641df7..f7e018104f9 100644 --- a/src/main/java/org/apache/sysds/common/Builtins.java +++ b/src/main/java/org/apache/sysds/common/Builtins.java @@ -96,6 +96,7 @@ public enum Builtins { CONV2D_BACKWARD_DATA("conv2d_backward_data", false), COOCCURRENCEMATRIX("cooccurrenceMatrix", true), COR("cor", true), + CORESETDT("coresetDT", true, ReturnType.MULTI_RETURN), CORRECTTYPOS("correctTypos", true), CORRECTTYPOSAPPLY("correctTyposApply", true), COS("cos", false), diff --git a/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinCoresetDTTest.java b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinCoresetDTTest.java new file mode 100644 index 00000000000..178ebaa30e9 --- /dev/null +++ b/src/test/java/org/apache/sysds/test/functions/builtin/part1/BuiltinCoresetDTTest.java @@ -0,0 +1,70 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +package org.apache.sysds.test.functions.builtin.part1; + +import java.util.HashMap; + +import org.apache.sysds.common.Types.ExecType; +import org.apache.sysds.runtime.matrix.data.MatrixValue.CellIndex; +import org.apache.sysds.test.AutomatedTestBase; +import org.apache.sysds.test.TestConfiguration; +import org.junit.Assert; +import org.junit.Test; + +public class BuiltinCoresetDTTest extends AutomatedTestBase { + + private final static String TEST_NAME = "coresetDT"; + private final static String TEST_DIR = "functions/builtin/"; + private static final String TEST_CLASS_DIR = TEST_DIR + BuiltinCoresetDTTest.class.getSimpleName() + "/"; + + private final static String WINE_DATA = DATASET_DIR + "wine/winequality-red-white.csv"; + private final static String WINE_TFSPEC = DATASET_DIR + "wine/tfspec.json"; + + // Test accuracy % drop of the coreset trained model vs the full data model + private static final double LOGREG_ACC_DROP_TOLE = 2.0; + private static final double DTREE_ACC_DROP_TOLE = 7.0; + + @Override + public void setUp() { + addTestConfiguration(TEST_NAME, + new TestConfiguration(TEST_CLASS_DIR, TEST_NAME, new String[]{"acc"})); + } + + @Test public void testWine10CP() { runCoresetDT(0.10); } + @Test public void testWine20CP() { runCoresetDT(0.20); } + @Test public void testWine50CP() { runCoresetDT(0.50); } + + private void runCoresetDT(double fraction) { + setExecMode(ExecType.CP); + loadTestConfiguration(getTestConfiguration(TEST_NAME)); + + fullDMLScriptName = SCRIPT_DIR + TEST_DIR + TEST_NAME + ".dml"; + programArgs = new String[]{"-args", WINE_DATA, WINE_TFSPEC, Double.toString(fraction), output("acc")}; + runTest(true, false, null, -1); + + HashMap acc = readDMLMatrixFromOutputDir("acc"); + check("logreg", acc.get(new CellIndex(1,1)), acc.get(new CellIndex(1,2)), LOGREG_ACC_DROP_TOLE); + check("dtree", acc.get(new CellIndex(2,1)), acc.get(new CellIndex(2,2)), DTREE_ACC_DROP_TOLE); + } + + private void check(String name, double full, double core, double tol) { + Assert.assertTrue(name + " coreset accuracy " + core + "% dropped more than " + tol + "pp below full accuracy " + full + "%", core >= full - tol); + } +} diff --git a/src/test/scripts/functions/builtin/coresetDT.dml b/src/test/scripts/functions/builtin/coresetDT.dml new file mode 100644 index 00000000000..2b1566d02b4 --- /dev/null +++ b/src/test/scripts/functions/builtin/coresetDT.dml @@ -0,0 +1,79 @@ +#------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +#------------------------------------------------------------- + +F = read($1, data_type="frame", format="csv", header=FALSE); +tfspec = read($2, data_type="scalar", value_type="string"); + +# as needed for decision tree +[D, meta] = transformencode(target=F, spec=tfspec); + +X = D[, 1:ncol(D)-1]; +y = D[, ncol(D)]; + +[XTrain, XTest, yTrain, yTest] = split(X=X, Y=y, f=0.8, cont=FALSE, seed=42); +ctypes = cbind(matrix(1, 1, ncol(X)), matrix(2, 1, 1)); + +[Xc, yc] = coresetDT(X=XTrain, y=yTrain, fraction=$3, seed=42, verbose=TRUE); + +t0 = time(); +Bf = multiLogReg(X=XTrain, Y=yTrain, verbose=FALSE); +tLrFull = (time() - t0) / 1e9; + +t0 = time(); +Bc = multiLogReg(X=Xc, Y=yc, verbose=FALSE); +tLrCore = (time() - t0) / 1e9; + +[Mf, pf, accLrFull] = multiLogRegPredict(X=XTest, B=Bf, Y=yTest); +[Mc, pc, accLrCore] = multiLogRegPredict(X=XTest, B=Bc, Y=yTest); + +# decision tree: fit on full vs coreset, evaluate on the test set +t0 = time(); +Tf = decisionTree(X=XTrain, y=yTrain, ctypes=ctypes, max_depth=12, min_leaf=5, min_split=10, max_features=1.0, seed=42); +tDtFull = (time() - t0) / 1e9; + +t0 = time(); +Tc = decisionTree(X=Xc, y=yc, ctypes=ctypes, max_depth=12, min_leaf=5, min_split=10, max_features=1.0, seed=42); +tDtCore = (time()-t0)/1e9; + +accDtFull = mean(decisionTreePredict(X=XTest, ctypes=ctypes, M=Tf) == yTest) * 100; +accDtCore = mean(decisionTreePredict(X=XTest, ctypes=ctypes, M=Tc) == yTest) * 100; + +realisedFrac = round(1000 * nrow(Xc) / nrow(XTrain)) / 1000; +tFullLr = round(100 * tLrFull) / 100; +dtimeLr = round(100 * (tLrCore - tLrFull)) / 100; +accFullLr = round(100 * accLrFull) / 100; +daccLr = round(100 * (accLrCore - accLrFull)) / 100; +tFullDt = round(100 * tDtFull) / 100; +dtimeDt = round(100 * (tDtCore - tDtFull)) / 100; +accFullDt = round(100 * accDtFull) / 100; +daccDt = round(100 * (accDtCore - accDtFull)) / 100; + +print("coresetDT-test: realised fraction = " + realisedFrac + "%"); +print("coresetDT-test: logreg full = " + tFullLr + "s | delta Core time = " + dtimeLr + "s | full acc = " + accFullLr + "% | delta acc = " + daccLr + "%"); +print("coresetDT-test: dtree full = " + tFullDt + "s | delta Core time = " + dtimeDt + "s | full acc = " + accFullDt + "% | delta acc = " + daccDt + "%"); + +acc = matrix(0, 2, 2); +acc[1,1] = accLrFull; +acc[1,2] = accLrCore; + +acc[2,1] = accDtFull; +acc[2,2] = accDtCore; +write(acc, $4);