This section describes the benchmark dataset used to evaluate the performance of the proposed model. The skin lesion images are part of the dataset, and they are downloadable from the International Skin Imaging Collaboration (ISIC) repository. The dataset contains two classes: malignant and benign tumors, along with respective image annotations provided by medical doctors who performed biopsies to diagnose them. Dermoscopy imaging was performed to better visualize tumor structures. The dataset contains 3,297 images, of which 1,800 are benign and 1,497 are malignant. The dataset is publicly available on Kaggle. Figure 2 displays some sample images from the malignant class, and Figure 3 displays some sample images from the benign class. The images are scaled down to 224 × 224 pixels.

The Stem Block serves as the foundational feature extractor, as shown in Figure 4. It processes input images into a lower-resolution, high-channel representation suitable for subsequent layers (Sandler et al., 2018). The Stem Block consists of an initial 3 × 3 convolution with a stride of 2, followed by batch normalization and a ReLU6 activation, forming the conv1_relu layer (Zhu et al., 2024).

Depthwise separable convolution is a computationally efficient alternative to regular convolution and is utilized extensively in lightweight neural network frameworks. As opposed to conventional convolution where the same convolution filter is applied for all input channels, it performs the operation in two distinct steps: depthwise convolution and pointwise convolution. This considerably alleviates the cost of computation and improves model efficiency without significantly degrading performance. In the depthwise convolution process, one filter is applied to each input channel independently, instead of summing up all channels together using several filters. Equation 1 represents the formula for the computational cost of depthwise convolution, where h is the height of the input feature map, w is the width of the input feature map, c is the number of input channels, and k is the size of the convolutional kernel.

Equation 2 represents the computational cost of pointwise convolution. In pointwise convolution, 1 × 1 convolutions are employed to fuse the features of all channels, effectively combining channel information. Here, d is the number of output channels.

Equation 3 represents the overall computational cost of depthwise separable convolution, which is the sum of the depthwise and pointwise convolution costs.

One of the key reasons for choosing the Stem Block is its complexity vs. performance trade-off. The architecture is tailored to improve feature representation through linear bottlenecks and expansion layers such that it can learn rich spatial and channel-wise information efficiently. While it is lightweight, it has been demonstrated to produce excellent classification accuracy on large-scale databases such as ImageNet, where it outperforms the majority of conventional architectures based on an efficiency-to-performance ratio.

Equation 4 represents the structure of an inverted residual block, where the input x is initially subjected to a 1 × 1 convolution to increase the number of channels, followed by a 3 × 3 depthwise convolution, then compressed using another 1 × 1 convolution. Linear activation occurs before residual connection addition. This design assists in maintaining the input while facilitating economical feature transformation at lower computational expense. In addition, it has also proven to be highly suitable for medical image classification tasks, such as the ISIC skin cancer dataset employed in this study (Surya et al., 2023).

The expansion layer is a key component in the architecture of deep learning models. Its purpose is to expand the number of channels before applying depthwise convolution. This operation allows the network to learn richer and more complex feature representations. The expansion layer increases the number of channels by a factor t, which is set to 6. The input channels C_in are multiplied by this factor to produce a higher-dimensional intermediate representation known as the hidden channels C_hidden. Equation 5 represents the number of channels after expansion. It captures finer details and richer spatial features.

Medical imaging demands models that effectively extract relevant patterns at a computationally reasonable cost, particularly with big datasets. Its low-cost feature extraction mechanism (Sawal and Kathait, 2023) makes it a very good fit for transfer learning. This guarantees that the model retains richer feature representations along with fitting itself to the sensitivities of medical image classification (Tika Adilah and Kristiyanti, 2023). In comparison with heavyweight architectures, it provides equivalent classification accuracy but with much reduced computational requirements, which makes it a perfect fit for large-scale processing. In addition, its strong ability to generalize well across different tasks and datasets aligns with the scope of this study, where accurate and efficient classification of images of skin cancer is of high importance (Ragab et al., 2022). Due to the advantage of its lightweight design, powerful feature extraction ability, high efficiency, and excellent performance in medical image processing, it is a perfect model for this research (Ekmekyapar and Taşci, 2023). Its computation cost at the expense of meaningful feature extraction makes it perfect for skin lesion classification, finally improving diagnostic accuracy and clinical decision-making (Xu and Mohammadi, 2024).

In this study, we used t-SNE to identify the most learnable feature representations from the model's intermediate layers. t-SNE is a non-linear dimensionality reduction method that maintains local structures in high-dimensional data by probabilistic similarity matching (Van der Maaten and Hinton, 2008). t-SNE first computes pairwise similarities between high-dimensional feature vectors. For two features f_i and f_j from a layer l, it calculates how similar they are using a Gaussian kernel. Equation 6 expresses the similarity of point j to point i as a conditional probability, where f_i and f_j are feature vectors in high-dimensional space, σ is the bandwidth of the Gaussian kernel, and ∥fi-fj∥2 is the squared Euclidean distance between the features.

To make the similarity symmetric and ensure consistency between i and j, pij=  pi|j+ pj|i2N is defined, where N is the number of samples. t-SNE maps high-dimensional data into a low-dimensional space using a student's t-distribution, which is a heavy-tailed distribution that helps prevent crowding, with one degree of freedom. Equation 7 represents the similarity between the points y_i and y_j, which are low-dimensional points, and the denominator ensures that all q_ij sum to 1 to make it a proper probability distribution.

As a next step, t-SNE minimizes the difference between the high-dimensional and low-dimensional similarity distributions using the Kullback–Leibler (KL) divergence. Equation 8 is the objective function to minimize the difference, where p_ij is the similarity between the points i and j in the high-dimensional space, and q_ij is the similarity between the points i and j in the low-dimensional embedding. This cost function encourages similar points in the high-dimensional space to stay close together in the low-dimensional map. It allows us to visually interpret clusters, patterns, and outliers in complex datasets.

Given the complexity of skin lesion classification, effective feature selection plays an important role in improving classification performance (Yuan et al., 2024; Villa-Pulgarin et al., 2022). Instead of arbitrarily choosing an intermediate layer for feature extraction, we utilized t-SNE to analyze the clustering behavior of feature representations from multiple layers within the pre-trained model's early activation block. Feature maps were extracted from several convolutional and dense layers. Each layer's output was reshaped into a feature vector representation. The extracted high-dimensional features were projected into a two-dimensional space using t-SNE. The objective was to observe how well the feature representations of benign and malignant skin lesions were clustered. Layers exhibiting well-separated clusters of benign and malignant lesions were considered the most learnable, as they provided high discriminative power. We selected the Conv1_relu layer from the Stem Block as it produced a well-clustered plot, as shown in Figure 12. The features from the layers were selected for further processing in the meta-learning classifier.

The Conv1_ReLU layer in the Stem Block consists of a convolutional operation followed by batch normalization and a ReLU6 activation function. The first layer of the network is responsible for extracting initial low-level features from the input image. It starts with a standard 3 × 3 convolution that processes the input RGB image X ∈ R^{224 X 224 X 3}. The convolution uses a filter tensor W ∈ R^{3 X 3 X 3 X 32}, which consists of 32 filters, each of size 3 × 3 × 3, along with a bias term B ∈ R³². Equation 9 defines the operation as follows:

where ^* denotes the convolution operation with a stride of 2. This operation reduces the spatial resolution and expands the depth channels from 3 to 32. Following convolution, the output is passed through batch normalization to stabilize training and improve convergence. For each output channel K ∈ {1, ... ,32}, the feature map is normalized using the channel-wise mean μ^(k) and standard deviation σ^(k), as shown in Equation 10. Equation 11 denotes that the normalized value is scaled and shifted using learnable parameters γ^(k), β^(k).

After normalization, the layer applies the ReLU6 activation function, which clips the values between 0 and 6. It is defined as ReLU6(y) = min(max(y,0), 6) and ReLU6(y) = min(max(y,0), 6). This helps maintain stability in low-precision computations. Equation 12 shows the final output of this layer.

Meta-learning is realized as a stacking classifier, which is a form of ensemble-based meta-learning. We do not perform task-level learning (similar to few-shot learning), but instead model-level learning, where the meta-learner is trained to optimally combine the predictions of many base classifiers. Each base learner captures a different facet of the input feature space extracted from the pre-trained deep model; their outputs are given to a higher-level learner (the meta-classifier), which learns to produce the final prediction. This also turns the model into some sort of a stacked architecture, enabling the model to utilize diverse classifiers with complementary strengths to enhance model robustness and generalization capabilities. This approach is common in standard meta-learning methods, as it uses knowledge across a range of models to improve classification performance, especially when noise or unbalanced data are present, such as in our classification of skin lesions.

Meta-learning is a method in which predictions of different base models are combined with a meta-classifier (Itoo and Garg, 2022), as shown in Figure 5. The meta-classifier is trained to weigh the base model outputs optimally and combine them for better predictive accuracy.

Equation 13 shows how the meta-classifier computes the prediction from all the base models. Here, B is the number of base models, y^b is the prediction from the base models, and f-meta is the final estimator, which is a Lasso-regularized logistic regression in this case.

In a meta-learning classifier, prediction is performed through a two-tier process involving multiple base models and a meta-model. The base layer consists of multiple individual models f_b, each trained independently. For every input X, each base model generates its own prediction y^b, as shown in Equation 14. Here, b indexes each base learner, and f_b(X) denotes the prediction made by the b-th model on input X.

The predictions from all base learners are then passed to a meta-layer. It learns to combine them optimally. This is performed using logistic regression as the meta-model, with L1 regularization. Equation 15 shows the final prediction being class 1, given all base predictions. Here, y^b is the prediction from the b-th base model, Bw_b is the weight assigned to each base model's output by the meta-model, w₀ is the bias term in the logistic regression, and σ() is the sigmoid activation function used to produce the final probability. To prevent overfitting and increase sparsity, an L1 regularization constraint is placed on the weight vector w, such that ||w||₁ ≤ λ. Here, ||w||₁ represents the L1 norm of the weights, and λ is a hyperparameter controlling the strength of regularization.

In the meta-classifier, six base models are employed—namely, CatBoost, XGBoost, SVM, random forest, gradient boosting, and k-nearest neighbors (KNN). CatBoost is a gradient boosting decision tree that is specifically used for categorical data. It is used for effective and precise classification and regression. It encodes the categorical features with target statistics from previous rows (Prokhorenkova et al., 2018). Equation 16 is the encoded value for the categorical feature of the i-th sample. Here, y_j represents the true target values of previous rows, and i-1 ensures that only preceding samples are used, preserving causality and avoiding target leakage. Equation 17 represents the objective function that includes both the loss term and regularization term. Here, n is the total number of samples, L (y_i, F(x_i)) is the loss function, F(x) = ∑t=1Tft (x) is the model's cumulative prediction over T decision trees, and Ω(f_t) is a regularization function applied to each individual tree f_t.

XGBoost is another gradient boosting algorithm that minimizes the loss function with a second-order Taylor expansion (Raihan et al., 2023). It is specifically built for high-speed classification and regression tasks, tuning both speed and accuracy. Equation 18 represents the objective function at the t-th iteration. Here, n is the number of training samples, L(yi, F(t)(xi)) is the loss between the true label and prediction at iteration t, T denotes the number of leaves in the newly added tree, and γ and λ are regularization hyperparameters that penalize complexity. Equation 19 represents the regularization component that encourages simpler trees. Here, γT discourages the number of leaves in the tree, discouraging overly complex models, and λ2∥w∥2 is a regularization term applied to the leaf weights w, controlling their magnitude and enhancing generalization.

SVC identifies a hyperplane that maximizes the margin between two classes. SVC supports linear and non-linear classification with kernel functions. It performs well in high-dimensional spaces and is resistant to overfitting, particularly in situations with a clear margin of separation between classes. Equation 20 represents the dual optimization objective performed by SVC. Here, α_i are the Lagrange multipliers to be improved, and K(x_i, x_j) is the kernel function that maps the data into a higher-dimensional space to make it linearly separable. K(x,z) = e^{−γ||x − z||2} is the radial basis function (RBF) kernel, where γ controls the width of the Gaussian kernel and determines the influence of each training example. This kernelized formulation allows SVC to handle complex, non-linear decision boundaries efficiently while maintaining strong generalization performance.

A random forest classifier is a machine learning model that employs many decision trees to classify data, making its predictions based on the majority vote of the individual trees. It is known for its accuracy and stability. It performs this by selecting a random subset of features and training the decision trees on a different subset of the data. The final prediction is made through majority voting (Pal, 2005). Equation 21 represents the prediction rule. Here, ŷ_RF is the final prediction from the random forest, ŷ_treek is the prediction from the k-th decision tree, and mode denotes the majority voting over all trees k = 1…K. Equation 22 represents the sampling correlation formula that is used to find how much the trees in the forest correlate in their predictions. Here, T(x; Z) is the output of a decision tree trained on the data subset Z with random feature sampling θ. Var_z and E_z denote the variance and expectation over the dataset subsets.

The gradient boosting classifier is an ensemble learning algorithm that constructs multiple decision trees sequentially, where each new tree corrects the errors of the previous ones. It reduces the loss function by optimizing residual errors using gradient descent. Gradient boosting combines weak learners to create a strong predictive model. Equation 23 denotes the prediction function at iteration m. Here, F_m(x) is the updated model after the m-th iteration, h_m(x) is the weak learner added at stage mmm, which fits the negative gradient of the loss function, and γ_m is the step size or weight for the m-th learner.

The KNN classifier is a machine learning algorithm that is non-parametric. It classifies an instance based on the majority class among its k closest neighbors in the feature space. The algorithm depends on distance metrics to compute similarity between points (Cunningham and Delany, 2021). Equation 25 represents the prediction rule of KNN. Here, N_k(x) denotes the set of KNN of instance x, and y_j represents the label of the neighbor x_j. Equation 25 provides the formula to calculate the distance between the points. Here, x_ix_j are two points in a p-dimensional space.

A meta-learner using logistic regression with Lasso regularization is an advanced ensemble learning technique, where logistic regression serves as the final classifier in the meta-learning model. Lasso (L1) regularization is applied to enforce sparsity by shrinking less important coefficients to zero, effectively performing feature selection. Equation 26 represents the objective function. Here, p_i = σ(w₀ + w₁x₁ + ... + w_px_p), Lasso penalty (||w||₁||w||₁) encourages sparsity, ŷ_i is the prediction from base learners, σ is the sigmoid function, w is the weight vector, ||w||₁ is the L1-norm of the weights, and λ is the regularization strength.